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ABSTRACT 

An analytical model is formulated for a three-dimensional nonlinear 
stability problem in a rocket motor combustion chamber. The chamber is 
modeled as a right circular cylinder with a short (multiorifice) nozzle, 
and an acoustic liner covering an arbitrary portion of the cylindrical 
periphery. The combustion is concentrated at the injector and the gas 
flow field is characterized by a mean Mach number. The unsteady combus- 
tion processes are formulated using the Crocco time lag model. The re- 
sulting equations are solved using a Green r s function method combined 
with numerical evaluation techniques. The influence of acoustic liners 
on the nonlinear waveforms is predicted. Nonlinear stability limits and 
regions where triggering is possible are also predicted for both lined 
and unlined combustors in terms of the combustion parameters. 
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constant defined in Equation 14 
speed of sound 

constant defined after Equations 18-2, 24 
constant defined after Equation 18-2 
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mass flow rate 

interaction index, longitudinal wave number, unit 
outward normal 

time dependent pressure 

time independent pressure 
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chamber temperature 
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chamber perturbation velocity 
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INTRODUCTION 

Combustion instability is generally classified into three types* Low 
frequency instability (ten to three hundred hertz) is usually characterized 
as a coupling between the propellant feed system and combustion chamber. 
High frequency instability (seven hundred to ten thousand hertz) is the 
most destructive and is usually characterized by an in-phase coupling be- 
tween the combustion processes and the pressure oscillations in the com- 
bustion chamber. The frequencies of these oscillations are close to the 
acoustic frequencies for the combustion chamber. Intermediate frequency 
instability, thought to be generated by a coupling of entropy waves and 
pressure oscillations, fills the frequency gap left by the two previous 
types of instability and is not frequently observed. This investigation is 
directed toward a better understanding of high frequency combustion in- 
stability. 

A rocket motor can be broken down into four areas of interest, the 
combustion processes, the nozzle, a damping device, and the gas dynamic 
field in the combustion chamber. The combustion processes, even in steady 
operation, are extremely complex and are not well understood. The pro- 
cesses involve injection, atomization, mixing, vaporization, and chemical 
kinetics. A hueristic approach to modeling the combustion processes in a 
gross way useful to the engineer was suggested by Theodore von Karman in 
1941i His basic hypothesis was that the combustion processes could be 
collected to form a time lag which represented the residence time of the 
propellants from injection to combustion. This early work led to the 

Crocco sensitive time lag theory which formulates a relationship between 
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the thermodynamic variables of the chamber and the time lag for an element 
of propellant. The time lag, T, represents the residence time before com- 
bustion for which the propellant element is exposed to gasdynamic oscilla- 
tions. The effect of the gasdynamic oscillations is correlated by the 
interaction index, n. This parameter is an indicator of the propellant 
elements T sensitivity to the chamber oscillations / 2,3,4 ' 1 This theory has 
been experimentally verified and the results show good agreement with the 
theory for the cases studied/^ 

The combustion parameters (n and t) are used to determine stability 
limits for a combustion chamber. Neutral stability limits separate the 
stability map into two regions. The first of these is one in which oscilla- 
tions grow in time and are termed unstable, in the second one the oscilla- 
tions decay in time and are stable. A combustor with combustion parameters 
that lie on the neutral stability curve are termed neutrally stable. 

The behavior of oscillatory flow through the convergent portion of a 

nozzle has received considerable attention. Tsien^ 6 ^ treated the linear 

small amplitude, problem for a one-dimensional, isothermal nozzle flow. 

( 4 ) 

Crocco extended the analysis by dropping the isothermal assumption. 

Crocco and Sirignano^ ^ later extended the analyses to consider three- 
dimensional flow in nozzles of varying geometry. The problems of finite 
amplitude oscillations (the nonlinear problem) has also been analyzed for 
three-dimensional oscillations of moderate amplitude by Zinn and Crocco/ 9,10 ^ 

The use of a damping device is of considerable importance in rocket 
engines. These devices may consist of baffles, acoustic resonator damp- 
ing devices, or a combination of the two types* Baffles are blades that 
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are placed perpendicular to the injector face. The mechanisms by which 
baffles work is not understood, and they must be treated as part of the 
gas dynamic problem within the combustor. The baffle problem will not 
be included in this report. Acoustic resonator damping devices are a 
category of devices which may be attached or machined into the cylindri- 
cal periphery of the combustion chamber. The resonant chamber for these 
devices may be modeled as Helmholtz resonator, quarter-wave, or half-wave 
tube. The theoretical analysis for small amplitude oscillations in 
Helmholtz resonators was done by Lord Rayleigh The theory has been 
extended to include nonlinear influences and empirically corrected through 
the years, * 9 Recent analyses have theoretically evaluated differ- 

ent resonator devices to include such effects as finite amplitude oscilla- 
tions, mean and oscillatory chamber flows, liner-mean- thro ugh flow, flow 
in liner backing distances (for Helmholtz geometries), and differences in 

the mean gas properties of the combustion chamber and the resonanting 

j . (16,17,18,19) ^ 

device. Experimental results show good agreement with the 

above theories, 

The gas dynamics of the flow in the combustion chamber is dependent 

on the combustion processes, the nozzle, and the acoustic absorber (and 

the baffles). By choosing particular models to describe the combustion, 

nozzle, and acoustic resonator behaviors, the stability of the flow in 

the combustion chamber can, in many cases, be determined. 

The stability of small amplitude oscillations in combustion chambers 

has been extensively investigated using the time lag model to represent 

(4) 

the combustion process. Early studies concentrated on the influence 
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ofthe combustion process and nozzle interaction on the gas dynamics. 

( 21 ) 

Oberg and Kuluva investigated the influence of a full length liner 
on a combustion chamber with respect to frequency and decay rate. How- 
ever, their work did not include the influences of mean flow in the cham- 
ber, of combustion, or of the nozzle. A later full length liner investi- 

( 22 ) 

gation was made by Priem and Rice which included the influences of 

mean flow, combustion concentrated at the injector and the nozzle. A 

similar invetigation to the one above with the combustion distributed 

( 23 ) 

axially was performed by Sirignano. These investigations showed that 

acoustic liners change both the stability of wave oscillations and their 
frequency of oscillation. These studies were fundamentally sound, but 
they were not realistic. In practice an acoustic liner does not cover 
the entire surface of a combustor but only a fraction of it. The linear 
analysis of the partially lined combustion chamber for both concentrated 
and distributed combustion, with mean flow, liner and nozzle influences 
has recently been completed at Colorado State University . This 
latter analysis determined the influence of partial length acoustic liners 
at arbitrary locations and of arbitrary length of neutral stability limits 
for the combustor. 

The major failure of the linear analyses just described is their 
limitation to very small amplitude oscillations and their inability to 
predict waveforms similar to experimental results. Allied with this is 
their inability to predict the possibility of triggered instability. 
Triggering is defined as the initiation of finite amplitude oscillations 
by the introduction of a disturbance of sufficient magnitude. This may 



- 5 - 


occur accidentally, through the combustion process, or intentionally by 

* 

the introduction of a high pressure gas pulse or an explosive charge. 

While it has been shown that linear analysis produces stability limit 
results that are in agreement with experiment, the results do not accu- 
rately predict the waveforms. The experimental waveforms are found to 
have sharp peaks and shallow valleys, or even discontinuities which the 
linear analyses do not model at all. 

Due to the complex nature of the problem the investigation of the 
nonlinear (finite amplitude) instability is not as fully developed as the 
investigation of linear instability. The behavior of finite amplitude 

transverse waves in a circular cylinder was investigated by Maslen and 
/ o a ^ 

Moore. Their investigation included only fluid mechanical effects 

and did not include the influences of combustion, mean flow, nozzle, or 
damping devices. They concluded from their work that transverse waves 
do not steepen into shock waves as longitudinal waves must* Their results 
are of interest to the field of combustion instability because the pre- 
dicted nonlinear waveforms are similar to those obtained experimentally in 

/ 

rocket engines. 

The success of the sensitive time lag concept in linear analyses has 

led to its extension into the nonlinear regime. The nonlinear analyses 

may be broken down into two groups. The first group studied the behavior 

( 97 ) 

of one-dimensional longitudinal oscillations. Sirignano J assumed con- 
centrated combustion and that the chamber was terminated by a ,f short M 
(quasi-steady) nozzle. He demonstrated the existance of periodic, finite 
amplitude, longitudinal waves near the linear stability limit. These 
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solutions were unstable with respect to the linear stability limits 
indicating the possibility of triggering. This analysis was extended 
by Mitchell^ to include the possibility of discontinuous wave forms. 

He included both concentrated and distributed combustion and, again, he 
used the "short" nozzle approximation. Finite amplitude, periodic, dis- 
continuous oscillations were shown to exist for both types of combustion. 

It was also demonstrated that intrinsically stable waves may be triggered 

( 29 ) 

to produce unstable waveforms, A. later analysis by Zinn and Lores v J 
with distributed combustion and "short" nozzle was developed to include 
the transient behavior of a longitudinal oscillation. The authors were 
able to show the transient and limit cycle behavior and waveforms, and 
found that in some regions the waveforms exhibit shock wave characteristics. 
They were also able to show that the final form of the oscillations was in- 
dependent of the initial disturbance. 

The logical extension of the one-dimensional longitudinal analysis 
to three dimensions was made by Zinn. His model assumed concentrated 

combustion and was terminated by his previously discussed three-dimen- 
sional nozzle model. Zinn showed £he existence of three-dimensional, 
finite amplitude, periodic oscillations. He was able to prove the possi- 
bility of triggering unstable oscillations, however due to the compli- 
cated nature of the final form of his equations, he was unable to determine 
tr i&£ er ing amplitudes. Powell using a different method of analysis 
studied essentially the same problem. For both concentrated and distri- 
buted combustion with the chamber terminated by a short nozzle he was 
able to obtain solutions to his equations and to determine nonlinear 
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stability limits. One main drawback in the two previous analyses is that 
they are forced by their respective analytic techniques to assume a solu- 
tion in one of the spatial directions. In other words, they assume a 

priori the final waveform for one of the dimensions, A more complete 

survey of the work in linear and nonlinear combustion instability may be 

(23 31 32 33 ) 

found in any of the excellent survey papers, * * 9 

From the list of previous accomplishments in the area of nonlinear 
instability, it is evident that it is desirable to find an analytical 

solution to the problem in which no approximations need be made about the 

form of the spacial waveforms and to include the influence of partial 
length acoustic liners on the final waveforms. This is, indeed, the goal 
of this analysis. 

In the present work the combustion process will be represented by 
the Crocco time lag model. In addition the combustion will be taken to 
be concentrated at the injector surface. This latter assumption is 
equivalent to saying that the region in which combustion takes place is 
small when compared to the length of the chamber, which has been shown 
experimentally to be relatively accurate in some cases. The model 

will include mean flow effects on the oscillations and the chamber will 
be terminated by a "short" nozzle. Since it is not the purpose of this 
analysis to design tuned acoustic liners for a given chamber configura- 
tion the acoustic liner boundary conditions will be represented simply, 
through an admittance coefficient. The use of the admittance coefficient 
is consistent with the previous work in the field since results of theory 
and experiments on the effect of acoustic liners are typically given in 
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terras of this parameter. It is an outgrowth of the early work in acous- 
tics and its extension to the area of combustion instability. The axial 
location and the length of liner will remain arbitrary in the analysis. 

The Green’s function technique used herein has been successfully 
applied to acoustic problems many times in the literature (for example 
Ref. 35). One of the first demonstrations of its applicability to the 
combustion instability was made by Culick v ; in his study of solid pro- 
pellant rockets* The solution technique used herein closely parallels 

(21) (24) (25) 

that used by Oberg, Mitchell , et al. and Baer, et al. This 

technique is extended to determine an analytical solution for the non- 
linear equations with discontinuous boundary conditions. 

In summary, this report presents a formulation of the three-dimen- 
sional nonlinear stability problem (subject to the aforementioned model 
assumptions), solves the resulting equations using applied analysis and 
numerical techniques, predicts the influence of acoustic liners on the 
nonlinear waveforms, determines nonlinear stability limits in terras of 
the combustion parameters n and T, and predicts regions of operation 
where triggering is possible for both lined and unlined combustors. 

THEORY 

In order to facilitate the modeling of the flow in a combustor cer- 
tain simplifications must be made in order to obtain governing equations 
which can be solved. The combustor model (Fig. 1) chosen for this analy- 
sis has the combustion zone concentrated at the injector, a uniform cylin- 
drical cross section, a finite Mach number mean flow in the combustor. 
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and is terminated by a "short" nozzle. An acoustic liner of arbitrary 
location and impedance may cover a portion of the cylindrical periphery, 
while the remaining uncovered chamber surface is assumed to be an acous- 
tically hard wall, that is to have an infinite impedance. 

The concentrated combustion assumption implies that the combustor 
is long compared to the region where the combustion takes place, or that 
the combustion zone is immediately adjacent to the injector surface. 

This assumption serves to separate the gasdynamics from the combustion 
processes and thus leads to considerable analytical simplification in 
the gasdynamic governing equations (discussed later on) because of the 
fact that the combustion processes appear only as a boundary condition 
and not in the governing equations themselves. The response of the com- 
bustion zone mass generation rate according to the Crocco sensitive time 
(4) 

lag (n-x) model, is assumed to be pressure dependent only. 

The short nozzle or multi-orifice nozzle^ approximation (a col- 
lection of small nozzles in the combustor exit plane) assumes that the 
individual nozzle dimensions are small compared to the chamber dimensions 
and the through-flow time is small compared to the period of the chamber 
oscillations. Therefore at each instant the individual nozzle behaves 
as if the flow were steady or that a constant Mach number condition 
exists at the entrance to the nozzle. 

The damping influence of acoustic liners for non-vanishing oscilla- 
tion amplitudes is a result of the formation of jets at the exit of the 
passage connecting the acoustic liner f s resonant cavity and the combus- 
tion chamber and the eventual dissipation of the jets' kinetic energy in 
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the combustor and the cavity, and to the frictional effects in the connect- 
ing passage. The effectiveness of an acoustic liner is a function of fre- 
quency, and maximum damping is obtained around the resonant frequency of 
the liner where the pressure amplitude in the resonant cavity is at a 
maximum. The resonant frequency of the acoustic liner is a function of 
the state variables and the physical dimensions of the acoustic liner. In 
this analysis the acoustic liners will be assumed to be tuned for all fre- 
quencies. That is, it will be assumed that the physical dimensions of the 
damping device will be adjusted so that the resonant frequency of the liner 
corresponds to the frequency of the combustor or that the impedance is 
independent of the chamber frequency. 

The flow of gas in the combustor is three-dimensional and it will be 
assumed the oscillations in the flow field are periodic and that one fre- 
quency is dominant in the chamber. The concentrated combustion model 
permits the gasdynamic field to be treated as a single constituent, source 
free gas. It is further assumed that the gasdynamic field is homentropic, 
irrotational , and calorically perfect. Transport phenomena, such as dif- 
fusion, viscosity, and heat conduction, are neglected in the governing 
equations. The absence of diffusion is consistent with the single con- 
stituent gas assumption. The absence of heat conduction and viscosity 
lead to conservative results since they result in damping of pressure 
oscillations. Two final forms of natural damping are neglected. Droplet 
drag is neglected because of the concentrated combustion assumption, and 
wall friction is neglected, which is consistent with the inviscid assump- 


tion. 
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Using the model described above, a mathematical representation is 
obtained by writing the conservation equations and the boundary equations 
for the combustor. 

The conservation equations are written as follows 
Continuity : 

3p* 5 + 

a~* + 7 • = o 


Momentum: 


p* 


DV* 

Dt* 


+ VP* 


0 


Homentropic Condition : 
p* = CP* Y 

State (perfect gas relation) : 
p* = p*R*x* 

The following scheme is now used to put the governing equations in 
a nondimens ional form. The state variables are nondimensionalized by 
their respective steady state or mean values. The characteristic velocity 
and time will be the mean of sonic velocity (a*) and the wave travel time 
(a*/R*) , respectively, while chamber dimensions are nondimensionalized 
with respect to the chamber radius. The dimensionless quantities are 


r* 

R* * 

Z* 

2 R* » 

t*R* 

^ 9 

a* 

p = 

P* 

P* 

y 

T* 

T = — , 

A 

&ll 

II 

> 


P* 

T* 

a* 
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The nondimens ional governing equations are 
Continuity ; 

+ V • (pV) =0 (1) 


Momentum: 


P 



Homen tropic Condition : 


P 



( 2 ) 


(3) 


State ; 

P » pT 


(4) 


The irrotationali ty assumption permits the representation of the 
velocity vector in terms of a potential function . 

V = V$ (5) 


By substitution of Equations (3) and (5) into Equations (1) and (2), it 
is possible to combine the conservation equations into a single partial 
differential equation in terms of the velocity potential. 


qv 2 $ - 


3 2 $ _ 3 
at 7 3t 


(vt 


7$) + (V$ • 7l>)V 2 <i> 


+ ~ V<& • • V$) + (Y-l)-ff- V 2 $ 


( 6 ) 


where Q is a constant of integration. The pressure is expressed as a 
function of the velocity potential as follows 


3$ .L 
3t 2 


Hi 

(Vl> * Vt) + p Y = Q 


( 7 ) 
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Lacking a direct solution of the above equations it is assumed that 
P and $ can be expanded in an infinite series in powers of some ampli- 
tude parameter, e, and that each order of the function $ may be repre- 
sented by a Fourier series in time, t. By making use of the periodicity 
and mean flow Mach number assumptions the velocity potential is then 
represented by 


00 j 

$ = Mz + Z Z 
j=l k=l 


eV j) 

9 (k) 


i (ku)> t 
e 


( 8 ) 


where M is the mean flow Mach number, z is the nondimens ional axial 
coordinate, i is the unit complex (/-T) , j is the ordering parameter, 

W is the complex chamber frequency, k is the frequency ordering param- 
eter, t is the nondimensional time, and <J> is the time dependent poten- 
tial function. 

The pressure is defined similarly 


P = 1 + Z Z p 5?? g 1 (kw)t 


j=l k=l 


(k) 


( 9 ) 


The chamber frequency is also represented by an infinite series in 
the amplitude parameter 


.=u/ 0) + Z eV*> 
j=l 


( 10 ) 


Before going further it is convenient to introduce a new variable 
y, where y is the product of the frequency and time 

y = cot (11) 

The boundary condition on Equation (6) is • n = £P where n is 


the outward normal for the appropriate surface, and $ is an average 
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specific acoustic admittance for that surface. The surfaces are defined 
as the injector plane, the nozzle exit plane, and the cylindrical periph- 
ery. The cylindrical periphery can consist of two or more surfaces , 
depending on the location of the acoustic liner. The acoustic admittance 
changes discontinuously with the change in surfaces. 

The injector boundary condition is derived using the time-lag model 
to represent the combustion processes. This model relates the rate of 
mass generation due to combustion to the steady-state mass flow as follows 


dTv 

m = m(l - 

dx ■ 

where — represents the rate of change of the time lag, and m is the 
at 

mass flow rate. This relationship assumes that the injection rate is con- 
stant at all times. The evaluation of will not be repeated here, 

because it is available and clearly presented in Crocco’s work.^ The 
final form of Crocco's results will be accepted as a postulate in this 
report. His relationship is 


dx 

dt 


n(P’(t) - P’(t - T)) 


where P 1 = P - 1. By arithmetic manipulation we obtain the following 
relationship between the admittance and the combustion parameters n and T. 


3 


i 



(P’(t) - P'(t— t)K 

P r (t) ’ 


( 12 ) 


Terms of p 12 and higher are ignored in this derivation. Nonlinear 

( 27 ) 

extensions to the above model have been formulated by Sirignano and 
Zinn. The linear formulation is used here because of its simplicity 



and because corrections to the primary combustion parameters n and T 
can be calculated through in a simple way. 

The chamber exit plane or nozzle entrance plane boundary condition 
is derived using the short or quasi-steady nozzle approximation. This 
approximation assumes that the dimensions of the nozzle are small com- 
pared to chamber dimensions. Such a quasi-steady condition can be approx- 
imated by a uniform distribution of many small individual nozzles over the 
chamber exit plane. The quasi-steady nozzle condition ^ is 




1 2y ; 


A{1 + 




k 2y ; 


JC t! . 

2(y-i) 


( 13 ) 


where A is a constant. 

On the cylindrical periphery the acoustic admittance need not be uni- 
form along the axial length, and when a partial length liner is present 
the admittance must change discontinuously. When a liner is not present 
or on the segment of the cylindrical surface which is not covered by a 
liner the acoustic admittance is assumed to be identically zero. This 
implies that the surface is an acoustically hard wall, that does not de- 
flect due to the pressure oscillations. The acoustic impedance is a 
coefficient which describes the interaction of sound waves with a locally 
reacting surface. The acoustic impedance for the absorber is related to 
the acoustic admittance by the following equation 




( 14 ) 


where k is the specific acoustic impedance for the absorber. The 
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acoustic impedance is the term generally used in acoustic literature for 
theory and experiment. In experiment (both acoustic and rocket test 
firings) it is easier to determine the impedance. The reacting surface 
in this case is the acoustic absorber. The impedance is a function of 
absorber dimensions and the frequency of oscillation within the absorber. 
There are many theories, linear and nonlinear, which can be used to cal- 
culate the impedance. However, the nonlinear and some of the linear 
theories use experimental data to define some of the parameters used in 
the theory, because of this a typical acoustic impedance is arbitrarily 
chosen and used. This impedance will be constant, which implies that as 
the frequency changes the absorber dimensions change such that the main 
chamber frequency and the absorber frequency are the same. The absorber 
dimensions are not of interest in this report, the reader is referred to 
the work of Mitchell, et al. for the procedure used in the designing 

an acoustic absorber for a combustion chamber. 

The expansions for p and oj in terms of e and harmonic con- 
tent are now substituted into the governing equations, Equations (6) and 
(7) and the result is separated according to powers of £ . In this work 
terms through 0(e 3 ) were considered. Similarly, the boundary condi- 
tions on the bounding surfaces are also expanded. Equations and appro- 
priate boundary conditions for 0(e), 0(e 2 ) and Q(e^) are thus obtained. 
These equations and boundary conditions are lengthy. They are written 
down completely in Appendix A. 

For the j ^ order in e there are j linear equations and boundary 
conditions for j harmonic components of 4 /^. That is, to 0(e) there 
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is a single linear equation for <J>^. To 0(e 2 ) there are two linear 


equations, one for and one for there are three 

(3) (3) (3) 

linear equations for ^(2) and ^(3)* 8 enera ^ f° rm of these 


, ( 2 ) 


equations is 


(j) 


7 2 4< J > - w"” 2 - F »»>) + F i«- 2) ) 

(k) U 3y 2 V*(k) J + ^N.L^ * 9 **•' 

where = cj>^ ^ * F L (9 (k) } is a Unear function of $ (kj and 

^Cj 2)) is a nonlinear function of lower order solutions. 


In order to convert the governing equations and their boundary condi- 
tions into integral equations , Green's functions corresponding to the 
appropriate frequency values are introduced. The Green’s functions satisfy 
the equation 


V 2 G (k) + (kw (0) ) 2 G (k) = <5(r/r o ) (15) 

where vS°° • n = 0 on all boundaries and 6 is the Dirac delta func- 
tion. The Green’s function is represented as an orthonormal eigenfunction 
expansion for the chamber with no mean flow, no combustion, and solid 
walls. The form of this expansion is given in Appendix A. 

The partial differential equations for all orders are then converted 

to integral equations using the appropriate Green’s function and follow- 

(34) 

ing standard techniques. Details of this transformation as well as 

the complete integral equations are given in Appendix A. 

To lowest order the linear integral equation is exactly the same as 
the one solved previously by Mitchell et al. Solution is obtained 

by successive approximation and results finally in a combustion zone 
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eigenvalue, (or and x^) and a solution for <j)^ of the 

following form 

<P (1) = $ + cos m0 2 Z y„ J (A. r) cos 
T Y . H £n m Jim L 

x n 


where $ is a solution obtained by separation of variables for the un- 
lined chamber, and is a two dimensional matrix determined by the 

successive approximation technique. 

To second order (in e) two integral equations result for the two 

( 2 ) ( 2 ) 

harmonic components <f> ^ and ^(2)' Solution of the integral equation 

( 2 ) ( 2 ) 

for 4* ( 2 ) determines the form of and the first order corrections 

to 3^ and to, and <o^. This solution is found to be (see 

Appendix A) <f> = 0* 0, to^ = 0. That is, to second order no 

change in frequency or combustion response eigenvalues is predicted, and 
the first harmonic component of the second order solution is identically 


The solution of the integral equation for results in the 

following expression for 


^(2) J Z 1 a Jlmn ^Jlm 
Z m n 


The three dimensional matrix is evaluated directly by integration 

of nonlinear functions of over the appropriate boundary surfaces 

and bounded volume of the chamber* No iteration or successive approxi- 
mation is necessary since <j>^^ is now a known function* Most of the 
integrals can be evaluated analytically, however the ones involving 
triple products of Bessel functions were evaluated numerically following 
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a technique discussed in Appendix B. 

To third order three integral equations for <f) , c})^ an d ^(3) 

may be written down. In this work only the first harmonic component 

( 3 ) 

for c(>^ was considered. The reason for this is that it is this equa- 
tion which serves to determine the second order corrections to the eigen- 

( 2 ) ( 2 ) 

values and co, and to . It is the determination of this 

eigenvalue relationship (which specifies the nonlinear stability behavior 
of the combustor) which is the main goal of this research. The third 
order corrections to the nonlinear waveforms are not determined, mainly 
from a profound awe of the computational complexities involved. Thus, 
the third order equations are invoked only to determine the nonlinear 
eigenvalue corrections; the waveforms are determined only to second order 
accuracy. The determining third order eigenvalue relationship is found 
to be (see Appendix A) 

</AK°v i> + Mif (1> ] d » 1 

a i 

= Z {0) ff% H(0 (1) , *[2)’ w (2) )ds Q -fff% F( 4 > (1) , $[1], w (2) )dV o 

*8 V 

o o 

( 16 ) 

where both H and F are known nonlinear function of lower order solutions. 

( 2 ) ( 2 ) 

The method used to determine $ and uj is to specify a 
particular combustor configuration, to solve the linear equations for the 
waveform and the acoustic admittance at the injector, to solve the second 
order equation for the nonlinear waveform, and finally to solve the first 
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harmonic third order equation for the second order correction to the 
acoustic admittance at the injector. A demonstration of the success of 
the technique in predicting nonlinear waveforms and nonlinear stability 
limits in combustors with and without partial length acoustic liners is 
presented in the next section. 

RESULTS AND CONCLUSIONS 

The analytical method developed in this investigation can be used 
to predict pressure waveforms similar to experimentally observed wave- 
forms and to predict regions where triggering is possible for both un- 
lined and lined combustors. A particular combustor configuration will 
be used to demonstrate the effectiveness of the method. The combustor 
configuration chosen to demonstrate the technique is one in which several 
of the parameters chosen have adverse effects on the calculations. The 
combustor chosen has a long chamber (L/R = 2.7), a high mean flow Mach 
number (0.33), and an acoustic liner with strong damping influence (k = 5) 
covering the chamber walls from the injector to one— third of the chamber 
length. Many combustors are shorter, have a lower mean flow Mach number, 
and although the admittance coefficient for the liner may be larger, the 
liner typically covers a much smaller portion of the periphery. All 
these changes cause positive influences on the calculational rapidity 
and convergence of the technique. The configuration chosen thus serves 
as an upper limit, as far as severity of calculational problems are con- 
cerned, for a wide range of combustors of interest. 

After the model parameters are fixed the linear equation is solved. 
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A solution to the unlined chamber is found by Bernoulli's method (separa- 
tion of variables) and then the effect of the acoustic liner on the in- 
jector acoustic admittance and the waveform is determined by a modified 
Green's function technique using Picard’s method (successive functional 
approximations) . 

Once the linear injector admittance and waveform are determined the 

second order correction to the waveform can be determined. Application 

of the Green’s function method leads directly to a solution involving 

only the eigenfunctions (£0 and a coefficient matrix (B.,) . The coef f i- 

cient matrix (B^) involves integrals of the products of the linear 

solution (<J>^^) and of the eigenfunctions. To solve for the second 

order coefficients the integrals involving the spacial variables 0 and 

z are evaluated analytically. There is no closed form solution for the 

majority of the integrals involving the Bessel function and these are 

evaluated numerically (see Appendix B) . 

Finally, with the second order coefficient matrix calculated, the 

( 2 ) 

second order correction to the injector admittance (B^ ) can be deter- 

mined through the use of the eigenvalue relationship. Equation (16). Use 

of the expression (3^, = -M[n(l - e^ WT ) - then results in two equations 

( 2 ) ( 2 ) ( 2 ) 

relating the three quantities n , T and to . 

( 2 ) ( 2 ) 

In the calculations performed n and T were required to lie 

on the normal to the first order stability curve, for each point on that 

curve (in other words for each value of For a given value of £, 

( 2 ) ( 2 ) ( 2 ) 

then, the corresponding values of to , n , T and D, the normal 
displacement from the linear neutral stability curve could be calculated. 
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Results of such calculations for the combustor chosen with no liner 
are shown in Figure 2. The linear stability curve separates the region 
of stable oscillations (outside the curve) from the region of unstable 
oscillations (inside the curve). The stability curve is parameterized 
in frequency with frequency increasing from right to left. For conven- 
ience the ratio of the frequency to the first transverse acoustic fre- 
quency is used to define the frequency points along the linear curve. 

The linear stability curve has a low minimum value for n^ which im- 
plies that the combustor at this point is very sensitive to disturbances. 
The linear curve also exhibits mixed mode behavior in the range of 
interest for the combustion parameters. That is, there is a transition 
from first transverse at low frequencies to the combined first transverse- 
first longitudinal at higher frequencies. In the figure this transition 
appears as the loop with the region to the right of the loop being pre- 
dominately first transverse acoustic and the region to the left of the 
loop being the combined first transverse-first longitudinal acoustic mode. 

The various types of nonlinear stability behavior predicted will now 
be discussed, again using the frequency ratio as parameter. For frequency 
ratios less than about .93 all normal displacements from the linear sta- 
bility curve are found to be negative (inside the neutral stability curve) . 
A region of stable finite amplitude oscillations is consequently indicated 
inside the neutral stability curve for this part of the stability map. 

Small amplitude oscillations grow in amplitude until they reach the ampli- 
tude predicted for the n and T of interest. 

For frequency ratios between about .93 and .965 along the neutral 
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8tability limit normal displacements are found to be positive (outside 
the neutral stability curve) and to depend on £ approximately quad- 
radically. This type of behavior is shown in Figure 3 for frequency 
ratios of .94 and .95. The region outside the neutral stability curve 
for this range of frequency ratios is a region where triggering is pre- 
dicted. Though this region is a linearly stable region, finite amplitude 
oscillations are predicted from the nonlinear analysis. The conclusion 
to be drawn is that, though the combustor is stable to small amplitude 
disturbances, perturbations of sufficient amplitude may excite finite 
amplitude periodic oscillations. For example, at a point displaced .075 
units from the point on the neutral stability curve corresponding to a 
frequency ratio of .95, a disturbance of amplitude ,175 is required to 
trigger the finite amplitude oscillations (see Figure 3) . 

Also shown in this region are dashed lines corresponding to the 
locus of points requiring the same amplitude disturbance to excite oscil- 
lations. The further away from the neutral stability curve one of these 
constant amplitude lines is located, the more sensitive is the location 
on the n,T plane to triggering. In particular at the frequency ratios 
of approximately .93 and .965 all lines of constant s pass through the 
neutral stability curve and no triggering is possible. 

Proceeding outward from the neutral stability curve at approximately 
a 45° angle from the point corresponding to a frequency ratio of .965 is 
a solid-dashed line. This line represents the approximate location of a 
triggering limit. Along this line disturbances of very large amplitude 
(beyond the amplitude for which this analysis may be considered accurate) 
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are necessary to trigger finite amplitude oscillations. Below this curve 
no triggering is predicted. 

For values of the frequency ratio of from .965 to 1.01 both positive 
and negative normal displacements are predicted. This behavior is seen 
in Figure 4. For a given positive normal displacement two, one, or zero 
finite amplitude solutions may be found. When two solutions are predicted 
the lower amplitude solution is assumed to be unstable, the upper stable. 
Then, disturbances with amplitudes less than that of the lower solution 
decay while those with amplitudes greater than that of the lower solution 
either grow (if they are below the upper solution) or decay (if they are 
above the upper solution) until this amplitude is the same as that of the 
stable (upper) solution. The lower solution serves as a limiting minimum 
amplitude for triggering in this case. The limiting displacement for 
this type of behavior occurs when the lower and upper solution coincide. 
The locus of these limiting displacements serves as a triggering limit in 
that, for displacements beyond this locus, no triggering is possible. 

The locus is shown as a solid-dashed line in Figure 2. Note that this 
line is very close to the neutral stability limit for this region. 

Negative displacements are also predicted in this region and are 
presumed to indicate the existence of stable finite amplitude oscilla- 
tions of relatively large amplitude (see Figure 4) . The dashed continua- 
tions of the various curves in Figures 3 and 4 indicate regions of 
sufficient amplitude that higher order corrections, not considered here 
may be quite important. 

The frequency regime from 1.01 on corresponds to the mixed mode 
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region mentioned earlier. Here again a triggering region along with 
lines of constant amplitude (dashed lines) and a triggering limit solid- 
dashed line are shown. 

The results for the combustor with an acoustic liner are presented 
in Figures 5, 6 and 7 and are similar to those of Figures 2, 3 and 4 
except that the stability curves are shifted to larger values of n. 

This indicates that a substantial overall stabilizing effect is produced 
by the liner. It can be seen in Figure 5 that the low frequency trigger- 
ing regime (frequency ratio less than about .95) is relatively smaller 
when a liner is present. In addition the lines of constant £ are more 
closely spaced Indicating that larger amplitude distrubances are necessary 
for triggering when a liner is present. This is also shown in Figure 6. 
(Compare with Figure 3.) 

The moderate frequency triggering region which was present without 
the liner is completely eliminated when the liner is in place. No two 
solution or two normal displacement region occurs. 

In the mixed mode high frequency region the size of the triggering 
region is reduced and the required triggering amplitudes are increased 
when a liner is present as was the case for the lower frequency trigger- 
ing region. 

The conclusion here is that the addition of an acoustic liner not 
only increases the linear stability of the chamber (the stability curves 
are shifted upward) but also decreases substantially the sensitivity of 
the combustor to triggering. The liner thus improves both the linear 
and nonlinear stability of the combustor. 
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The nonlinear waveforms were determined by adding the nonlinear 
waveform perturbation to the linear waveform. In Figure 8 a typical 
experimental pressure oscillation (from Ref. 23) is shown. The oscilla- 
tion displays a sharp peak with a distorted valley. Figures 9 and 10 
are the pressure waveforms as a function of time at the injector for two 
typical frequencies for an unlined chamber. In Figure 9 the combustor 
frequency is below the resonant frequency and in Figure 10 the combustor 
frequency is above the resonant frequency. Figure 9 shows some nonlinear 
behavior in the slight steepening of the wave with a slightly sharper 
peak and shallower valley. Figure 10 has a nonlinear waveform with 
steepening of the waveform at the peak and a shallow valley. The results 
in Figure 10 are similar to these presented in Maslen and Moore v in 
their acoustic work, and the results presented in Powell and Zinn^^ in 
their nonlinear stability work. 

For a lined combustor the pressure-time plots (Figs* 11 and 12) are 
similar to those for the unlined chamber. In Figure 11 the nonlinear 
waveform coincides with the linear waveform. The liner has a sufficiently 
strong influence at this value of frequency so that there is no distor- 
tion. At the higher frequency (Figure 12) the influence of the liner is 
not as strong and in fact the wave steepening is more pronounced than in 
the unlined chamber. In Figures 13-16 pressure is plotted versus the 
axial dimension. Figure 13 shows that at a frequency below the acoustic 
frequency for an unlined chamber the nonlinear wave shows steepening and 
a sharper peak than the linear wave with an increase in the maximum and 
a decrease in the minimum pressure amplitude. For the same chamber at a 
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frequency above the acoustic frequency (Figure 14) the linear wave is a 
very flat first longitudinal wave. The nonlinear wave shifts to a more 
pronounced waveform similar to a second longitudinal wave with a decrease 
in the maximum pressure amplitude. Figure 15 shows the nonlinear wave 
exhibiting typical nonlinear behavior with steepening, shallower valleys, 
and sharper peaks for a lined chamber with increases and decreases in 
the maximum and minimum pressure amplitudes. Figure 16 is the lined 
chamber at a frequency above the acoustic frequency. The nonlinear wave 
exhibits typical nonlinear behavior with a decrease in the maximum per- 
turbation pressure. 

The analytical technique presented in this report is quite general 
and from the results the following conclusions are drawn: 

1. A Green’s function technique can be successfully applied 

in the solution of nonlinear combustion stability problems, 

2. The technique can predict experimentally observed nonlinear 
phenomena . 

a. regions where finite amplitude disturbances may be 
trigger instability in a linearly stable engine. 

b. nonlinear pressure waveforms similar to those experi- 
mentally observed, 

3. The technique is applicable to a wide range of combustor 
configurations as determined by the ability to handle the 
"limiting" case used. 

4. The addition of an acoustic liner to a combustor improves 
both the linear stability (higher combustion zone response 
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is required) and the nonlinear stability. (Regions where 
triggering is possible are reduced and larger amplitudes 
are required.) 

5. The technique requires a reasonable expenditure of computer 
time to determine the nonlinear behavior of a given com- 


bustor behavior 
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Figure 2. Stability map for an unlined combustor. Chamber length =2.7, mean flow 
Mach number = 0.33, ratio of specific heats =1.2. 
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Normal Displacement 


Figure 3. Normal displacement for an unlined combustor 
for various frequencies along the linear 
neutral stability curve. 



Figure 4. Normal displacement for an unlined combustor 
for various frequencies along the linear 
neutral stability curve. 
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Figure 5. Stability map for a lined combustor. Chamber length = 2.7, mean flow 
Mach number = 0.33, ratio of specific heats = 1.2, X = 0.0, X- = 0.90, 
k = 5.0. 1 z 
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Figure 6. Normal" displacement for a lined combustor for 
various frequencies along the linear neutral 
stability curve. 
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Figure 7. Normal displacement for a lined combustor for 
various frequencies along the linear neutral 
stability curve. / 
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Figure 8. Experimental pressure waveforms versus time 
(at the injector) (see Ref. 23 , p. 188). 
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Figure 9. Pressure versus time for an unlined combustor 
below resonant frequency. 
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Figure 10. 


Pressure versus time for an unlined combustor above 
the resonant frequency. 





Figure 11. Pressure versus time for a lined chamber below the 
resonant frequency. 
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Figure 12. Pressure versus time for a lined combustor above the 
resonant frequency . 
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Figure 13, Pressure versus axial coordinate for an unlined combustor 
below the resonant frequency. 
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Pressure versus axial coordinate for an unlined combustor 
above the resonant frequency. 
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Figure 15. 


Pressure versus axial coordinate for a lined combustor 
below the resonant frequency. 
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Appendix A 

This Appendix presents the differential and integral equations, 
boundary conditions, and eigenvalue relationships developed for the 
application of this technique. 

Substitution of the expansions in £ discussed in the Theory sec- 
tion into the governing equations Equations (6) and (7) leads to the 
following equations to the appropriate orders in £ 
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oz ay (k) az 3y 

+ -|v 2 $ (1) (V% (1) • 7^ (1) )| + -| • f$ C1) ) 

+ 2 J°V« -^ + <a»«V« + J 1 ) 2 ) 2i|^l 

l2 A (2) /-, v 

+ 2 “ ll) + i a? <* a) • * a) >] + “ C2) « y^- 

<U> 


+ (y-l)[a) (1) V 2 $ (1) ^— ] 

»$ - ^ # > a)3 - ^ ♦ “ (0) ^ 
.(2) 


(3) 


(3) 


3$ 


+ .a> + „<*> ag. + *a> . *g>| 

The associated boundary conditions are written 


0(e): 

V& (1) • n - 3 (0) p (1) 

0(e 2 ) : 

v$( 2 > . i = b (0) p (2) 

v (k) n p p (k) 

0 (e 3 ) : 

rifcO) • n = ft(°)n (3 > 

V 0 (k) n " 3 P(k) 


At the liner 6 


( 0 ) 


B P( k ) + P p 

» 8^ = 6^ * 0, where k is the acoustic 


impedance and is taken to be a real constant. At the nozzle M, 
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3^ a 3 ^) a o, where M Is the Mach number. At the injector (com- 
bustion zone) 3^ , 3^\ and 3^^ are complex eigenvalues to be 
determined by the periodicity requirement. 

The partial differential equations along with their associated 
boundary conditions are converted to integral equations using a Green's 
function represented by an eigenfunction expansion as follows 


,00 _ 


W r)i WV 

r, ”“ n (to (0> )2-n? 


= s.r t . 


5.mn 


is an orthornormal eigenfunction given by 


fl 


_ cos m9 „ nrrz 

- ~M— J m ( \m r)cos ~ 
£mn 


and is the frequency of the closed chamber 


n? „ X 2 + . 

Jlmn lm t 2 


The quantity is a root of the equation 


J' (X 0 r) = 0 and A fl 
m Jim r=l Jlmn 


is the normalization constant which is determined from 


'"Am dV * 1 ‘ 

Utilization of this Green's function in the governing equation and 

(34) 

following standard techniques' ' converts the partial differential equa- 
tions to integral equations for each power of e and each harmonic com- 
ponent. To lowest order (0(e)) 
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( 1 ) 


4> (1) = // G^ 1) (r/r^)6[yiw (0) 4) (1) + yM d S 


8z 


+ m g* 11 (?/?„>[** + 21al «» M ]4 » o 

o 


( 1 ) 


(18-1) 


d) 


CO) 2 


»M - V* + 2i “ (0> M d ¥ c 

o 


( 1 ) 


3* 


+ // n N [ylw (0) (j) (1) + yM 3z 

s 

o 


(i) 


"] d S, 


(18-2) 


where Is one of the acoustic eigenfunctions for the chamber with no 

mean flow, no combustion zone, liner or nozzle, rv is the eigenvalue 
for the frequency of oscillation associated with the eigenfunction 
and G,^ is a Green’s function for the chamber with one term (d„) 

Vi N 

removed from the series for G ^ The solution to the linear equations is 
given formally by 


nlTz 


= 9 + cos 1 8 ijj E w £n J m^Jtm r ^ COS 2 


where 


1 cos m 0 T /, , izB . „ izB^ v 

4> = ~z — ■ — J~(x* £ r )( e + Ce 2 ) 

AK ^ ¥ m 1 m 
5, m n 

to (0) M ± [o) ( 0>2 M 2 + (M 2 -l)(Xg ~ - cn^ 2 )]* 5 


>1,2 


(1-M 2 ) 


- C B > + V<“ <0> + “l>l 

— — e A ■ ■ ■ ■■ y — tj- 

l b 2 + B n y(w to; + mb 2 ) j 


and ¥ is a normalization constant determined by 
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!! J * dV ' !• 

T 

The coefficients and the specific acoustic admittance at the in- 

jector are determined by a successive approximation technique (essen- 
tially Picard’s method). The subscripts f.,ra,n represent the radial, 
tangential, and longitudinal wave numbers, respectively. The reader is 
referred to the papers of Mitchell, et al.^ 2 ^ and Espander^ 37 ^ for a 


more detailed description of the above solution technique. 

To second order in e two harmonic components exist, <j> 


( 2 ) ( 2 ) 

4> ^ 2 ) * equation for is 


( 2 ) 

(1) 


and 


( 2 ) 


9 2 <J) 


( 2 ) 




N 


3z 


3z 2 


+ u U) (u ,»|^. u (0) (1) ](W 

o Z c 

a*< 2) 

r g™ <?/?„) { y b w) [m -5 

s 


+ SS g<«(?/? o h yS (0) [m + 


+ ^ * a) ] + YS a) [M 2|21 + iu W*0)]} ds . 


where the surface integral represents the integral over all three surfaces, 
This equation in turn is divided into an eigenvalue relationship for 
6^^ and and a linear homogeneous integral relationship for ^(2)* 

These are, respectively 


B<» - „<»{/// £» {0) * a> - dv 
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- ff Y - ^- - fl N (r o )fl (1) dS Q }/{// Yfl N (r o )(ito (0) ^ (1) 

s s 

o i 


+ « ^> ds i> 


aA U) a2(fl u; 

*(1) ■ S! J g n 1>( ' / ' o ,[21 “ <0> m -ap- + m2 

o 


+ ss ^’G/z o )yis w Cm -5^ + ioi (0 >*< 2 ’] ds o 
s 

o 


The only consistent solution to the eigenvalue relationship is 3^ } ■ 

“0. The solution to the integral equation for ^( 2 ) on t * ie ot ^ er 
hand is = a<}>^ where 0 is arbitrary. In particular a may be 

identically zero. Thus the first harmonic second order equations imply 

= 0) (1) - <j>^ = 0. 

(2) 

The integral relationship for ^ 2 ) 

*(2) = G i 2) [ f l«(2)) + 81 (^ (1) )]dV o + ff GjJ 2) [f 2 (cf,^> 


+ g 2 0t> (1) >]ds. 


fl(4>J?s) = 2i(2w (0) )M 4F + M 
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+ i (M ^- + i&J (0) ) (f<l> (1 ‘ > • vY l h 


-^[vV 1} • 4< l) - +M ^-) 2 ] 


3z 


The equation may be rewritten so that the terms involving a PP ear 

on the left hand side of the equation and the terms remaining are a 
function only of 


♦«) - ff J g n 2) £ '<*(y> d¥ o - " G N 2>f 2»S )ds o 


■ SSS G N 2, 81<+ U) > dv 0 + !S '4 2) S2(* (1) MS ( 

^ „ s 

o o 


The right hand side of this equation can be rewritten as R'.H.S, 

E 0 Z Z B„ Q„ , where 
X, m n Jlmn £mn 


^v.( r ) 

2 J /// - N ° 


^ 8lW<1))dV o 


+ rr fl N (r o ) ,,Uk q 

ff (0),2 82<(|> )dS o-l 


S o <(2u w ')‘ - Dp 


Assuming that a particular solution can be represented as a triple 

(2 ) ■> 

infinite series, <J>^ = Z^a^Cr) , substitution into Equation (20) yields 


Vj S k a ijk (1 + Y ijk + 6 ijk )n ijk (r) « E A E m S n B £*Amn <r) 


where 
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hjk - - °W ^iJk^O^O 


S! f 2 (S5ijk(r 0 ))dS o 

s 


Due to the orthonormality of the eigenfunctions, the coefficients 
a j£mn are ^ e ^^- nec ^ as 


Amn ‘ (1 + 5 *mn + W ; 

The integrals appearing in the coefficient B^, Y^* and 6^ can 
be evaluated analytically, except for the Bessel integrals which must be 
evaluated numerically. The Bessel integrals involve triple products of 
the Bessel functions of different orders. The evaluation of these inte- 
grals is discussed in Appendix B. 

To third order in £ three harmonic component equations occur. As 
discussed in the theory section, only the one for <f>^ is of interest 
since it yields the relationship between the eigenvalues and 

The equation for is 


♦8) ■ S! J «$ 1) <?/f o >r<a» <0) M 


+ FCr o )]dV o + /; G«> (t/t o )y 6 «»C( lu ‘ 0 )*<5> + M 


- e (0) H(? o ):ds o + ref' i G‘ 1) Cr/? o )[i 11 <0) *< l > + M-^idSj 



A-9 


where 

F(r) -i (Ho (0) +M^-)?* (2) • Vf (1> + (^ti) V 2 » (1) (M 

+ 2ic <0) * (2 >) + vV 2) ( m 4^-+ i«. (0) + (1) ) 

C/ Z 

+ j (v| (1) • V<F ( 1 ) )V 2 <f> (1J + -| (V$ (1) * V<f>^)V 2 <(> (1) 

+ -|vtj) (1) • . V$ (1) ) +-v| (1) • ^(V4> (1) • V<j> (1) ) 

+ „<«[» (0)^(1)] 

and 

H<;> - - $91 p(« 2 p«> - £*.«>♦<» + vi< 2 > • *«] 

Following essentially the same separation procedure as was used in 

the development of the second order eigenvalue relationship the following 

(2) (o') 

equation relating 6^ and or y results 

Y6< 2> //V^ (0) * (1) + » ^d Sj[ 

= g (0) //O. H(? )dS - Iff F (r )dV 

s N o' o N o' o 

o o 

where F(r^) and H(r Q ) are t ^ le functions defined previously. 



B-l 

Appendix B 

The integrals appearing in Equations 25 and 30 are Lommel-type inte- 
grals involving products of Bessel functions to the cubic and quartic 
powers. The integration scheme chosen to evaluate these integrals was an 
improved Gaussian Quadrature technique by Kronrod.^^ In Gaussian Quad- 
rature one gives up the freedom of arbitrary end-points (a,b) for the 
interval and fixes them at (-1,1) by a change of variables from x to y 
by x = -~(b-a) y + -|-(b+a) . Assuming the original integrand <J>(x) is a 
continuous function, a new integrand is introduced as a polynomial in y. 
b i 

/ <Ji (x)dx = S <j>(y)dy = <x <f> (y ) + aj^Cyj) + ... 

a -1 ° 0 

where $(y). = y^^Cy) and L ft (y) is the Legendre polynomial. The 
are evaluated (and corresponding values of y^) by evaluating <p(y) 
for increasing values of n and evaluating the weighting functions 
(cO and the nodes (y^) . The technique is accurate up to a polynomial 
of order exactly 2n - 1. Kronrod’s method is similar except that a 
second orthogonal function is introduced and the weights and nodes are 
evaluated. 

b +1 

I $(x)dx - f <J> n (y) ^ n+1 Cy> d y * 3 0 <Ky 0 ) 

a -1 

where “ Y Q . + Yi* 1 + y 2 X 2 + ... + This quadrature method 

has an accuracy of 3n+l for even n and 3n+2 for odd n. 

The forms of the integrals to be evaluated are 
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1 

/ 

0 




V X £m r) ^ a im r)dr 


( 1 ) 


1 

/ r n J (X e r) J ft (Aft ft r) J-(Ar-r) J , (A 0 , ,r)dr 
n a w m X,m m £m m' ' Jl'm’ ' 


( 2 ) 


where -1 < n (an integer) £ + 1 
0 m, &, m, m’ £ 2 


J’(A 0 r) 
m Jim ' 


- 0 


r=l 


The values for were obtained from tabulated values in Refer- 

ences 38, 39, and 40. 

Three studies were made to determine the accuracy of the Improved 
Quadrature and the number of nodes needed. Fettis^ X ^ evaluated the 
integral 


1 

/ c J l^ X JLm r ^ dr for P ar£icular values of the 

subscripts using Gaussian Quadrature. His results are published to five 
place accuracy with an upper bound on the error of 0(10*"^). For most 
of the cases run nine nodes were the maximum needed to reproduce the 
same solutions to the published five decimal places. The second study 
involved a comparison using Simpson’s rule and Improved Quadrature. The 
Simpson's rule algorithm (Conte had an internal error control which 
was set at one tenth of one percent. The results for both integrals (1) 
and (2) were identical within the error limit. However the Simpson's 
rule algorithm was slower than Improved Quadrature by a factor of ten. 
The third check on convergence was an internal check made by increasing 
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the number of nodes and calculating the error with respect to the largest 
node case run (21 nodes). (This corresponds to a 65th order polynomial.) 
Depending on the integral evaluated the calculation time varied from 3 
to 21 milliseconds for 21 nodes on a third generation CDC computer. 

Using this internal comparison it was found that 9 nodes gave an error 
which was less than 0.001 percent when compared with the 21 node solution. 
The 9 node Improved Quadrature technique was chosen for evaluating inte- 
grals (1) and (2) . 
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Appendix C 

This section describes the use of the computer program, COMBADM, to 
determine regions of triggering and nonlinear stability limits for a 
given set of input parameters which characterize the combustion chamber. 
The main program calls three main subroutines LINEAR, NONLINE and EIGEN. 
LINEAR calculates the linear coefficient matrix, y, and the injector 
admittance, When the calculation for the linear solution are com- 

pleted all variables that are not necessary later are destroyed. 

NONLINE calls the routines which calculate the second order coeffi- 
cient matrix. On the call to NONLINE a Hollerith variable is used to 
determine whether or not the linear cross product terms are calculated. 
If the hollerith variable equals CROSSbPROD the cross product terms are 
calculated. For any other hollerith variable they are not calculated. 
After the calculations are complete the program returns to the main pro- 
gram. 

EIGEN is now called by the main program. This is the subroutine 

which calculates the second order correction to the injector admittance 

and the second order correction to the frequency, w . It should be 

( 2 ) 

noted that a numerical root solving technique is used to find co x y and 
since the equation is transcendental care must be taken in assuring that 
the proper root is found. Several numerical root solving techniques 
were tried and it was found that the Secant Method was the most 
dependable. The subroutine SOLSION is the root finding technique. 

The input parameters are: 

1. LENGNO - the problem identifier 

2. F - the initial frequency ratio of the combustor 
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3. 

4. 

5. 


6 . 


7. 


8 . 


9. 


10 . 


11 . 


12 . 


13. 


14. 


15. 


16. 


17. 


18. 


ALENGTH - 
AMACH 
GAMMA 
G 

AK 

XI 

X2 

ERR0R1 - 

ERR0R2 - 

LHAT 

MHAT 

NHAT ' - 
LDEX 

NDEX 

MATRIX - 
FINC 


the length to radius ratio for the combustor 
mean flow Mach Number 
ratio of specific heats 

complex number specifying the acoustic impedance 
of the nozzle 

complex number specifying the acoustic impedance 
of the liner 

the beginning of the liner opening with respect 
to the injector 

the end of the liner opening with respect to 
the injector 

the maximum acceptable error in the real part of 
the combustion admittance 

the maximum acceptable error in the imaginary part 

of the combustion admittance 

radial acoustic mode number 

transverse acoustic mode number 

longitudinal acoustic mode number 

specifies the number of radial terms in the 

coefficient matrix 

specifies the number of longitudinal terms in the 
coefficient matrix 

identifier used to determine whether or not the 

coefficient matrix is printed out 

the increment size of the frequency ratio 
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19. FMAX - the frequency for which a combustion admittance 

is calculated 

20. EPSILON - initial wave amplitude 

21. LOWLIM - minimum acceptable value of w 2 

22. OMEGING - increment of u) 2 to start root search 

23. UPLIM - maximum acceptable value of oj 2 

Sample Calculation 

An example of the output from COMBADM is shown on the page immediately 
following the listing of the program. It is suggested that the output be 
used for comparison to insure correct operation of the program. In the 
output the linear neutral stability point for the unlined combustor is 
denoted as BETAWI and the linear neutral stability point for the lined 
combustor to the specified accuracy is the last value of BETAI(J) listed. 
The n and T values listed next are given for the lined combustor first 
and the unlined combustor next. 

The nonlinear (second order) coefficient matrix is printed out next 

for the unlined combustor. The n and T values above the line follow- 

( 2 ) 

ing the coefficient matrix are for oi =0. The numerical values in- 
side the lines are the second order corrections. These are followed by 
the corrected values for the frequency, the admittance, n, and t. The 
normal displacement is the distance between the linear stability point 
and the stability point for the given amplitude parameter. The quanti- 
ties G and G1 are the error checks for the root u> 2 . Following 
the results for the unlined combustor are the results for the lined 



combustor which are presented in the same form as above 


The form of the input for COMBADM is shown immpediately following 
the list of the program. The input is on five data cards as follows: 


Card 1 

Columns 

Variable 

Type 

1-10 

LENGNO 

Alphanumeric word 

11-20 

Re(E) 

Real number 

21-30 

Im(F) 

Real number 

31-40 

ALENGTH 

Real number 

41-50 

AMACH 

Real number 

51-60 

GAMMA. 

Real number 

61-70 

Re (G) 

Real number 

71-80 

Im(G) 

Real number 

Card 2 

Columns 

Variable 

Type 

1-10 

Re (AK) 

Real number 

11-20 

Im(AK) 

Real number 

21-30 

XI 

Real number 

31-40 

X2 

Real number 

41-50 

ERR0R1 

Real number 

51-60 

ERR0R2 

Real number 

61-63 

LHAT 

Integer 

64-66 

MHAT 

Integer 

67-69 

NHAT 

Integer 

70-72 

LDEX 

Integer 

73-75 

NDEX 

Integer 
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76-80 MATRIX* Alphanumeric word 

Card 3 


Columns 

Variable 

Type 

1-10 

Re(FINC) 

Real number 

11-20 

Im(FINC) 

Real number 

21-30 

Re(FMAX) 

Real number 

31-40 

Im(FMAX) 

Real number 

Card 4 



Columns 

Variable 

Type 

26-35 

EPSILON 

Real number 

Card 5 



Columns 

Variable 

Type 

1-10 

LOWLIM 

Real number 

11-20 

0MEGINC 

Real number 

21-30 

UPLIM 

Real number 


* If the coefficient matrix is to be printed out, then MATRIX is printed 
as YESBB, where B signifies a blank. 



C-6 


co<°> 

«<2) 


M 

Y 

^£.mn 

<S„ 

Jlmn 






A 




Amn 


B 

B 

C 

i 


1 

2 


3< 2) 

»i 

^c 

^£ran 


Nomenclature 

COMB ADM 

OMEGA 

0MEGA2 

OMEGAC 

AMACH 

GAMMA 

Function GAMM 
Function DELTA 
Function ETA 
Function BESPRIM 
Function BESSEL 
Function FNORM 
Function PSI 

BETA1 

BETA2 

C 

I 

BET AW 1 

BETAI or BETA1S 

BETAI2 

BETA1C 

BETAC 

BETAN 

MUSND and MUCROSS 
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MUSND 

£mn 


AMU 

n 


N or B 

T 


TAU 

e 


EPSILON 

* 
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PROGRAM COMRAOM ( TAPE1 » TAPE?, TAPER, T APE6 • OUTDUT = 10l > 

REAL LENGTH . IJ1 , TJ? , IJ3 , IJ4 
REAL N , NVEC , NSAvFM) , T AUSAvE (3) 

COMPLEX PR I , FNORM , AN 

COMPLEX E » G« AK, EINC* EMAX, 8ETASTR 

COMPLEX BETA], BETA?, BETAN, BETAC, BFTAI, RETAWT, OMEGA. C. Ak \ , 

] AK?* PSI?. AMU, BETA I S * Cl?, I 
COMPLEX D 1 » 0? » D? , DA , D5 » D6 

COM 4 ON / «LKa / BET A 1 , HPT A? , AMACH. GAMMA, OMEGA, ALENGTH, C« AK) 
1, AK2, M~T-iN. ^LMOA, BET A I » BETAC, BFTA-VI, LHAT. MHAT. NHAT. 
1PSI?*X1,X? 

COMMON / BLKB / AMU ( 3, 30 . ?> . BET A I S ( ?>, I, MATRIX 

COMMON / BL*C / LENGTH , DlOO) , D?<3 A > • D 3 < 3 n J , DA(30) . 06( 30 

5) , D6 ( 30 ) . Al(?) 

COMMON / BLOCK / Cl ? (30 , 30 ) * I SET ,LDEX , NDEX . F.RR -.)P ) , ER>- 0* 

CO' MON / NORM / N(3> , TAU(3) ♦ NVEC , TAUVEC , PRINT 
C *********** ********** &»«***»»«***** **** *ff«* ****** **«»» *** 
c *** 

C *** THE FIRST THREE OATA CAROS ARE THE general engine characteristics 
C »** THE FOURTH DATA CARO IS FOR T HF POST . r ION AND CASF IN PRFSS AND VFL 
C *** THE NEXT OATA CARO S^ECIF T ES THE INITIAL TIME , THE TIME INCREM , 

C AND The EPSILON 

C *** 

C *** 

C <»*» L ENG NO IS THE PROBLEM IDENTIFIER 
C *** F IS THE COMPLEX FREQUENCY RATIO 
C »** ALENGTH IS THE CHAMBER LENGTH TO RADIUS RATIO 
C *** AMACH IS THE CHAMBER MEAN FLOW MACH NUMBER 
C *** GAMMA IS THF RATIO OF THF SPECIFIC HEATS 
C *** G IS THE COMPLEX N077LE AOMITTANCE 
C AK IS THE COMPLEX LINER ADMITTANCE 

C XI IS THE INITIAL LINER LOCATION 

C *** X? IS THE LOCATION OF THF F^D OF THF LINER 

C *** ERROR! IS THE MAXIMUM ALLOWABLE ERROR IN THE PEAL PART OF Thf S OLN 

C *** ERROR2 IS Thf MAXIMUM ALLOWARLE ERROR IN THE IMaG PART OF THF SOLN 

C ««* LHAT IS THE PRIMARY RADIaL MODE 
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C *** MHAT IS THE PRIMARY TANGENTIAL MODE 
C NHAT IS THE PRIMARY LONGITUDINAL MODE 

C *** LDEX ANO NOEX SPECIFY THF SI7E OF THE COEFFICIENT MATRIX IN Thc l 
C »** AND N DIRCTIONS RESPECTIVELY 

C **« IF MATRIX IS YES THF COEFFICIENT MATPICIES ARE PRINTED OUT 
C FINC IS THE FREQUENCY RATIO INCREMENT 

C **» FMAX IS THE MAXIMUM VALUF WHICH THE FREQUENCY RATIO m A y TAKE 
C **# EPSILON IS THF AMPLITUDE PARAMETER 

C *** LOWLIM IS THE INITIAL! ZTTON FOR THE OMEGA (?) SEARCH 
C *** OMEGINC IS THE FREQUFNCY INCREMENT FOR OMEGA (?) 

C UPLTM IS THE MAXIMUM VALUE THAT OMEGA (?) MAY TA«E 

C *** 

WPITE<6»1000> 

1000 FORMAT (*PM NEW RTQBON PL 7 *> 

2 READ (St?flO>LFNGNO»F, ALENGTH t AMACHtGAMMA*G»AK,Xl .X?«FPPORl TERROR?. 
1LHAT»MHAT«NHAT*LDEX,NDEX.MATRIX»FINC.EMAX 

IF < EOF • S ) SO * 3 

3 BETAC = I./GAMMA/AK 

G = CMPLX f (GAMMA + 1 .)/P./r,AMMA»0.0 ) 

BETAN = AMACH * 1 G - 1. / GAMMA > 

JX = 30 
JY = JX - I 
PRINT = 10H 

IF < BETAC .FQ* CMPLX (O.PtO.O) .OP. X? .EQ. 0.0 > JY = 0 
CALL LINEAR ( LENGNO • F • G. * AK • FINC ♦ FMAX , JX ) 

CALL NTAU ( 1 ,/GAMMA-BETAlS <1 > /AMACH , OMEGA ) 

NSAVE(l) = N S TAUSAVE(l) = TAU 

CALL NTAU ( I./GAMMA-RETAWI/AMACH . OMEGA ) 

N ( 3 > = N S TAU<3> * TAU 

F = F ♦ FINC 
JX = 30 

CALL LINEAR ( LENGNO t F * G * AK ♦ FINC * FMAX , JX ) 

CALL NTAU ( l ./GAMMA-SETA IS ( 1 ) /AMACH , OMEGA ) 

NSAVE ( ?) = N . $ TAUSAVF<?) = TAU 

CALL NTAU < 1 ./GAMMA-BETAWI/AMACH , OMEGA ) 



C-10 


N(2> = N <E TAU (?) = T At) 

CALL wRyTF 
F = F + FINC 

4 CALL LINFAR ( LENGNO * F « G 9 AK , FINC • FMftx . JX ) 
CALL NTAU { | . /GAMMA-BET A I S (1 ) /AMACH , OMEGA ) 

NSAVE<3) = N S TAUSAVEO) = TAIJ 

CALL NTAU ( l./GAMMA-RETAWI/AMACH , OMEGA ) 

CALL WRYTF 

CALL NORMAL 

CALL REED 

MATRIX = 3HYES 

LENGTH = ALFNGTH 

HETdl = BFTAWI 

0 ETASTR = BETAC 

BETAC = CMPLX ( 0.0 « 0.0 ) 

CALL INI T ( CI2 ♦ AMU(I*t«2> • 01 * 3 * 30 > 

WRITE (6.600 ) 

A I ( 1 ) = 3. 14IS9?6S35ftQ79 
A1 ( 2 ) = Al ( 1 ) / 2. 

WP I TE (6.601) 

CALL CONST 

CALL NONLINE ( 1 CHL INEAR ) 

CALL EIGEN 
GO TO 6 
6 CONTINUE 

IF ( JY .EO. 0 .OR. JX ,FO. 30 ) GO TQ 5 
BETAI = BE T A I S ( ? > 

BETAC = BFTASTR 
WR I TE ( 6.602) 

CALL CONST 

CALL NONLINE ( 10 HCROSS > 

N ( 3 ) = NSAVE( 1 > % TACK 3) = T AUSAVE ( 1 ) 

NSAVEU) = N(?> * TAUSAVE(l) = TAIJ (?) 

N { 2) = NSAVE (?) % T AU ( 2) = T AUSAVE ( 2 ) 

NSAVE (2) = N(l) 5 TAUSAVEt?) = TALJ(l) 

N ( 1 ) = NSAVE ( 3 ) $ TAU(l) = TAJJSAVEd) 
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CALL NORMAL 
CALL EIGEN 
N ( 3 ) = NSAVE<1> 
NSAVEU) = N<2> 
N<2) - NS A VE ( 2) 
NSA VE (2) » N ( 1 ) 
1 F = F ♦ FTNC 


% TAU(3) = TAUSAVEU) 
* TAUSAVE(l) = TAU (2> 
S TAU(2) = T AUSAVE (2) 
5 TAUS AVE ( ? ) = T AIJ( I ) 


JX = 19 

IF ( F ,GT • FMAX ) 2 * e> 

5 N ( 3) = N(?> S TAU ( 3) = TAtJ<2) 

IF ( JX .FQ« 30 ) GO TO «;(> 

N (2 > = N(l) % TALK 2> = TAU(l) I GO TO 1 

SO CALL EXIT 

200 FORMAT (A10 t7Fl0.0/6F10, 0 .SI 3 »A5/AF 1 0.0) 

201 FORMAT (2A10.3F10. 0.2110) 

202 FORMAT!//) 

600 FORMAT <*1*135(1H=)/* *46 ( 1H=) .43X . 46 ( 1 H= ) /» *46(1H=)* THE FOLLOWI 
SNG ARE THE NONLINEAR RESULTS *46<1H=)/* *46 < 1H=) ,<.1X.*6 < 1H=> /* *1 
fSSUHEJ) 

601 FORMAT (*0*50 ( IHa) /* SECOND ORDER. NO LINER */* * 50 (lHv)> 

602 FORMAT (*0*50(1Ha)/* SECOND ORDER. W/ LINER */* *S0<lHv>> 

END 
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FUNCTION BtSPRIM (H, L) 

DIMENSION A < 1 

THESE ARE the ROOTS OF the DERIVATIVE of THE BESSEL FUNCTION Or OR 
C#**» SET EQUAL TO ZERO 


A { 1 « 

i> 


0.00000000 

A (2 « 

l) 

— 

3.83170597 

A { 3 . 

1) 

= 

7.01558667 

A (A ♦ 

l) 

s 

10.17346814 

A (5 V 

1) 

JT 

13.12369194 

A < 6 * 

1) 

= 

16.47063009 

A ( 7 * 

1) 


19,61585891 

A < 8 * 

1) 


2?, 76008439 

A (9 ♦ 

1) 

— 

29.90367209 

A ( I A ♦ 

I) 

S 

29.04682893 


C**»* THESE ARE THE ROOTS OF THE DERIVATIVE OF THE BESSEL FUNCTION OF OR 


SET 

EQUAL 

TO ZERO 

A ( 1 

♦ 2) 

X 

1.84118379 

A (2 

»2) 

as 

5.33144277 

A (3 

1 2) 

z 

8.53631637 

A (4 

»2) 

— 

11. 70600490 

A <5 

*2) 

- 

14.86398863 

A (6 

♦ 2) 

— 

18.01592786 

A ( 7 

*2) 

=r 

21.16436996 

A (8 

♦ 2) 


24.31132686 

A (9 

» 2 ) 

- 

27.45705097 

A < I 6 

♦ 2) 


30.60192297 


THESE ARE THF ROOTS OF THE DERIVATIVE OF THF BESSEL FUNCTION OF OP 


SET EQUAL TO ZERO 


A ( 1 

.3) 

3 

3.09423693 

A (2 

931 

s 

6.70613319 

A (3 

♦ 3 ) 

= 

9.9694678? 

A (4 

»3> 

= 

13.17037086 

A (5 

*3) 


16.34792232 

A (6 

♦ 31 

rr 

19.51291279 

A (7 

*3) 

s 

22.67198177 

A (8 

♦ 3 ) 

= 

25.82603714 
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A ( 9 *3) = 28. 07767277 
A (10*3) = 32.1273270? 

C**«* THESE ARE THE ROOTS or THE DERIVATIVE OF THE BESSEL FUNCTION Of 0 

C«**« SET EQUAL TO ZERO 

A ( 1 «4) = 4.20118894 
A ( 2 *.4) = 8.11523660 
A <3 *4) = 11.34592431 
A ( 4 *4) = 14.58584829 
A(5 • 4 I = 17.78874787 
A (6 »4) = go. 97247694 
A ( 7 »4) = 24.1&489743 

A ( 8 *4) = 27.31005793 
A (9 *4) s 30.47026881 
A(l6*4> s 33.62694018 

C**«* THESE ARE THE ROOTS OF THE DERIVATIVE OF THE BESSEL FUNCTION OF 0 


c«««« 

SET 

equal 

TO ZERO 


A ( 1 

♦ 8) 

S 

S. 31755313 


A (2 

*5) 

= 

9.28239629 


A ( 3 

♦ 5 ) 

= 

12.68190844 


A (4 

♦ 5) 

= 

15.96410704 


A ( 5 

♦ 5) 

— 

19. 19602880 


A (6 

♦ 5) 

= 

22.4010 3227 


A ( 7 

*5> 

— 

25.58975968 


A (8 

»5> 

= 

28.7678362? 


A ( 9 

.5) 

= 

31.93853934 


A { 1 9 

*5) 

= 

35.10391668 

1 

besprim 


AIL. M> 


RETURN 

ENTRY BESSEL 

C***« THESE ARE THE BESSEL NUMBERS OF ORDER ZERO FOR Th£ 7F*0S OF THF f 
C*«** FUNCTION 


All 

* 1) 

= 

1.00000000 

A (2 

*1) 


-0.4027588095 

A ( 3 

*1) 

= 

0.301128303 

A (4 

♦ 1) 

“ 

-0 • 249704877 

A (5 

*1) 


0.218359407 
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A (6 *1) = -0,19645371 
A < 7 *1) = 0.150053375 
A (8 *1) - -0.167184600 
A ( 9 .1) = 0. 156724985 
A d 0 ♦ 1 ) = -0.148011108 

<3**** THESE ARE tHE BESSEL NUMBERS OF OPOF D ONE FOR THE 7E w OS OF T^E a 
c***» FUNCTION 

A(1 *2) = 0.5818649368 
A ( 2 *2) = -0,3461258542 

A ( 3 *2) = 0.2732981131 
A ( 4 *2) = -0,233304416 

A ( 5 *2) = 0.207012651 
A { 6 *2) = -0.188017488 
A ( 7 ,2) = 0.17345905O 
A (8 ♦ 2 1 = -0.161838211 
A < 9 .2) = 0.152282069 
AdO*2> = -0.144242905 

C***» THESE ARE THE BESSEL NUMBERS OF ORDER TWO FOR THE 7FROS OF THE R 
£«««« function 

A ( 1 * 3 ) = 0.4864961885 
A (2 *31 = -0.3135283099 

A ( 3 *3) = 0.2547441235 
A (4 *3) = -0.220881581 

A ( 5 *3) = 0.197937434 
A ( 6 *3) = -0.181010000 
A ( 7 * 3 J - 0. 167835534 
A ( 8 ♦ 3 ) = -0.157195167 
A (9 ,3) = 0.148363778 
A { 1 6 * 3 ) = -0.140878333 

C#w** THESE ARE tHE BESSEL NUMRERS OF ORDER THREE FOR THE 7EROS OF THF R 

C«mh»* function 

A(1 *4) = 0.4343942763 
A { 2 *41 = -0.2911584413 
A {3 *4) = 0.240738175 
A {4 *4) = -0.210865204 

A (5 *4) = 0.190419022 
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4(6 ♦<*) = -0.I7504840S 
A ( 7 *4) = 0.168954965 
A<8 »4) = -0. 15310240Q 

4(9 »4) = 0. 144866574 
A ( 1 f) ♦ 4 ) = -0.137844513 

THESE 4RE THE BESSEL NUMBERS OF OROE^ FOUR FOR THE ?FROS OF ThF 

c***» function 

A(L »S> = 0.3996S14S48 
A ( 8 *5 ) = — 0 . 5743809949 
4(3 * 5 ) = 0. 859590468 
A ( 4 *5 ) = -0.508763849 
A ( 5 *5) = 0.184059896 
A (6 »5> = -0.169878516 
A (7 *5) = n • 1 58655378 
A (8 *5) = -0.149451156 
A (9 *5) = 0. 141714307 
A { 1 0 * 5 ) = -0.I3S0863?8 
GO TO 1 
END 
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SUBROUTINE 8ETAIJ ( J > 

COMPLEX AI50»Al5l.AI52»AT53«A IS5 , A I 5^ J, OUM l«PSl,AMU»8ETAIS,CO M FGA» 
i DUMMY . I « C1P.SUM1.SUM2.SUM3.SUMA , AN 
COMPLEX BETA 1 * RETA2. BETAN. BETAC. BETAI. PET AW 1 4 OMEGA. C. A<t. 
1AK2. PSI2. FNORM, EPSIR7M 

COMMON / BLKA / RET&1. B^TA2, AMACH, GAMMA. OMEGA. ALENGTH, C. A*1 
1. AK2 . RETAN. RL M OA • RET A I ♦ RETAC. «ETAWI. LhaT* MHAT • . NHAT . 
1PSI2.X1.X? 

COMMON / RLK0 / AMU ( 3t3"*2) .RETAIS ( 2). I. MATRIX 
COMMON / RLOC.K / C 12 ( 30 . 1ft.) , I SET . LDEX . NOEX . ERROR 1 .ERP0R2 
NHAI = NHAT ♦ 1 

C0ME6A=CMPLX {-A I MAG {OMEGA) .REAL (OMEGA ) ) 

SUMl=SUM2=SUM3*SI|MA=CMPLX (0.0.0. 0) 

OO 10 LP = 1 . LOEX 
OUM =BESSEL (MHATM .LP) 

DO 10 NP = 1 ♦ NOEX 

SUMa*SUM4 ♦ AMU(LP.NP» I ) *DUM *C 1 2 (NHAT* 1 .NP) 

10 CONTINUE 

00 20 NP * t ♦ NOEX 
DUM1=AMU(LHAT.NP» 1> 

SUM 1 =SUM 1 ♦ OUM 1 

SUMP=SUM2 ♦ DUMI*{-J ,)**(MP *1> 

IF (NP.EQ.NHAT +1) GOTO 20 

SUM?=SUM3 ♦ OUMl*2.*COMEGA*(NP-l>**2/HNP-l) # *2 -NHAT**2> 

1*(1. -(-l.)**(NP .NHAT =))) 

20 CONTINUE 

A150*HETAI*AKI/PSI (DUMMY) /EPSIZN (LHAT, MHAT .NHA T , ALENGTH , RLMOA > 

A IS ) =BETAN*GAMMA*COMFGA* { (-1 .) **NHAT) *SUM2*FN0RM (LH AT.MH A T ♦ NH AT . 

1 RLMOA. ALENGTH >/EPSl?N(LH AT. MH AT, NHAT .ALENGTH, RLMOA) 

A 152 - ( * 1 ♦ ) * A MACH* SUM 3*F NORM (LHAT .MHAT .NHAT ♦ PLMDA • ALENGTH) /EPSI7N( 
1LHAT.MHAT, NHAT, ALENGTH, RL.MOA) 

A 1 51=RETAC*FN0RM(LH AT.MH AT, NHAT, RLMOA, ALENGTH) *BFSSEL(MHAT + 1 .LHAT) 
i»SUMA-/EPSIRZN (LHAT .MHAT .NHAT .ALENGTH. RLMOA) 

AI5S= CAK1/PSI (DUMMY) + G uMMA»COMEGA*FNORM (LHAT . MHAT , NHAT .RLMOA . 

1 ALENGTH) »SUMl) /EPSI7N (LHAT, MHAT, NHAT. ALENGTH, RLMOA) 

BETAIS (2)= (AISO -AIS1 - AI52 - AIS3)/AIS5 
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t 
j 

WRI TE < 6 ♦ 1 00 ) J« BETAI5(2) 

AN = l* / GAMMA - BETAIS(2> / AMACH 
WRITE 101) AN 
A 150 = BETAIS (2) - BET A 15 ( 1 ) 

IF ( AHS (REAL (AT50I/REAL (RETATS ( 1 ) ) ) .LE. ERROR] . AND. ABS ( A IMAG ( A IR 
10)/AI M AG(RETAIS( I > ) ) .L^.ERROR?) GO TO ? 

HETAlSd) = BETA T S ( ?> 

rc-tkrn 

2 I SET = 1 

F = REAL ( OMEGA ) / l«8Atl0178 
RETURN 

100. FORMAT (*0 BETA I (*♦ I2«*) = *»?G21.1A) 

101 FORMAT (*♦*« 70K** N = * « 2G21 • l A » * I* > 

ENO 
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SUB»OUTINE CALMIM JP ) 

COMPLEX AMU.G1 .G2.PST .DUMMY. I ,DUM3,C1? 

COMPLEX OljM i .OUM2* BFTAIS » SUMS . S * COMEG A , SUM1 *SU M ? , SU M 3 » SUM4 
COMPLEX BETAl* BETA?* BETAN* BETAC. BETAl. BETA W I ♦ OMEGA* C* A^l* 
1AK2. PS 1 2* FMOPM 

COMMON / PLK A / BETAl. BFTA2, AMACM* GAMMA * OMEGA* ALENGTH, C* AKJ 
i* AK2* BETAN* RL M OA ♦ BET A I * RETAC, RET AW I * LHAT* mhaT* NHAT. 

1 PS1? * X 1 ♦ X? 

COMMON / BLKR / AMU ( 3 * 30 *2) . BET 0 IS < ?>. I, MATRIX 
COMMON / BLOCK / C 1? HO * BO) ♦ I SET *LDEX ,NOFX , ERROR 1 .EPBOR? 

ISET = 0 

DO BO L = 1 * LDEX 
RLMHA1= 8ESPRIM(MHAT*1*L) 

AC= BESSEL (MHAT»1 »LHAT) /RESSEL (MHAT*1 * L) 

00 SO NX = l « NOEX 
N=NX-1 

IF (L.EQ.LHAT.ANO.N.EO.NHAT) GOTO 40 
DUM\s CMPLX(0.*0*> 

IF (L *EQ.LHAT ) DUM1 = <BETAT-BETAWI ) * AK \/ ( EPS 1 7N (LH A T *MH4T *N« ALENGTH. 
iRLMDAM 

OUM 1 - (BETAC* (G1 (X?.N> -fi\ ( Xl *N) ) /EPSIK AT (LHAT*MHAT*N* ALENGTH. RLMi) 
IA 1 ) »AC«-0UM1 ) 

0UM?=OUMl/( (OMEGA**?-ETA(LHAT*MHAT.NX. ALENGTH, RLM0A1) ) *PSI ( DUMMY ) ) 

GOTO 50 

40 DUM? = CMPLX ( 0.0 *0.0 ) 

50 AMU (L.NXfl)=DUM2/ENOPM(LHAT*MHAT,NHAT*RLMOA* ALENGTH) 

COMEGA=CMPLX (-AIMAG(OMEGA) .REAL (OMEGA ) ) 

NX = l 

OUMl=CMPLX (0.0,0. 0> 

OUM?=CMPLX (0.0.0. 0) 

N = NHAT 

NX = NHAT ♦ 1 

DO ? NY = 1 . NDEX 

np=ny- 1 

IF (N.EO.NPi B.7 

7 C12(NX.NY)=GAMMA/2./(N*»P-NP*»2)*(ALENGTH*C0MEGA/3.14l5R2G tt ( (N*NP> 
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1* (SIN( <N-NP) *3. l&159?6*X?/ALENGTH)-5rN ( (N-NP) *3. 1415926*X 1/ALENGTH 
1 ) ) ♦ (N-NP)* (SIN( (N*NP>*3. 1415926*X2/ALENGTH) -SIN( (N*N») *3. 141S9?6*X 
U*ALENGTH) ) )-AMACH*N=>*( (^P»N)*<COS( (NP-N)*3. 14159?6*X2/ALENGTH)-cO 
IS ( (NP-N) *3. 1415926*Xl/ALFNGTH) ) ♦ (NP-N) *(COS( (NP*N> *3 . 1M5924*X?/AL 
1ENGTH)-C0S< <NP+N)*3,1415Q26*Xl/ALENGTHn )) 

GO TO ? 

a IF(NP*EO.O)R »10 

9 C12<NX»NY)s(X?-X1)*C0MEGA*GAMMA 
GO TO 2 

10 C12(NXfNY)sGAMMA*( ( ( S IN <P. *N*3 . 141 59?6*X 2/ALENGTH > -S I M ( 2. *N*3 . 1 4 1 5 
1926*XI/ ALENGTH) ) /A ./N/3. 141 5926* A LENGTH* (X2-X1 ) /?. ) *C0MEG A-AM AfH/p 
1.*<SIN (N*3, 1415926*X2/ALFNGTH) **?-SIN(N*3. 1 4 ] 59?6*X 1 /ALFNGTH ) **2 ) ) 

2 CONTINUE 

CALL dETA I J( 2 > 

IF(JP.EO.l) RETURN 
00 l NX s I , NOEX 
N = NX - 1 
00 1 NY = 1 . NOEX 
np=ny-i 

IF (N»EQ»MP) 3,4 

4 C12 <NX,NY)=GAMMA/2./ (N**?-NP*»2) * (ALFNGTH*C0MEGA/3« 14159P6* ( (N*NP> 
1*(SIN( (N-NP1*3.1415926*X?/ALENGTH)-SIN( (N-NP)*3. 1415P?6»X1/ALENGTH 
1 ) ) ♦ (N-NP) * (SIN ( (N*NP)*3. 1415926*X2/ALENGTH> -SIN ( (N*NP> *3. 14J5926*X 
l l*ALENGTH> > ) -AMACH*NP* ( (MP+N) * (COS ( (NP-N) *3, 1415926*X2/ALEN&TH) -CO 
iS< (NP-N) *3. 1^159?6*Xl/ALFNGTH> ) + (NP-N) * (COS ( (NP*IM> *3. 1415926*X?/AL 
lENGTH) -COS < (NP*N) *3» 14 15Q?6*X1 /ALENGTH) ) ) ) 

GO TO 1 

3 IF(NP,EQ. 0>S«6 

5 C12(NX,NY)=(X2-X1)*C0MEGA*GAMNA 
GO TO 1 

6 C12<NX,NY>=GAMMA*< ( ( SIN ( ?. *N*3 . 1 4 15926* X2/ALENGTH) -$ IN ( 2. *N*3. 1 4 15 
1926*X1/ALFNGTH) ) /4. /N/3. 1 4] 5926*ALENGTH+ ( X?-X 1 ) /? . ) *COMEGA-AMACH/? 
1.*(SIN (N*3. 14 15926*X?/ALFMGTH)**2-S IN (N*3. 14 159?G*X1 /ALENGTH) **2) ) 

1 CONTINUE 
DO 60 J=2,JP 
DO 30 LS s 1 , LOEX 
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00 «0 NX s 1 » NDEX 
NS=MX-1 

PLM0A2*RES p PlM ( MHAT ♦ 1* LS) ,.,. T 

QUH1 = 8ETAC*E p SIZN (LHAT.MHAT« NS , ALENGTH.PLMOA) *HES5EL MHAT* l . 
1LHAT) /EPSIK AT (LS.MH AT. NS. ALENGTH.PLMDA2) /BESSEL (MHAT* l.LS) 

OOM l = CMPLX ( 0.0* 0.0 ) 

SUMS=CMPIX (0. .0.) 

IFC LS .EQ. LHAT .AND. N<= .EO. NHAT ) 15. 15 

15 AMU(LHAT. NX . ? ) * CMPLX ( 0.0 * 0.0 ) 

GO TO 80 

16 SUMl«SUM2eSWM3*SUMA*CMPLX <0.0 *0.0) 

DO 90 LP = 1 . LOEX 
RLMDAT=8ESSEL (MHAT *l»LP) 

RLM0A4=BESSFL (MHAT *1,LhaT> 

DO 90 Ny = 1 . NDEX 

SOMA=SUMA * AMU(LP.NY. 1 ) *C 1 2 (NX . N V ) 

90 CONTINUE 

SUMA=SUM4*DUM1*RLMDAT/RLM0A4 
00 ?1 NP = 1 . NDEX 
DUMt=AMU<LS»NP» 1) 

SUM 1 =SUM 1 *DUMl 

SUM?=SUM2 ♦ DUM1*(-1.>**(NP *I> 

IF (NP . EQ .NX ) GOTO 21 

SUM3*SUM3 ♦ OUH1*2.*COMEGA*<NP-1)**2 *< 1 .-<-1 . > ** <NP*NX> ) /< 

1 (NP-1)**2 -NS**2 ) 

21 CONTINUE 

SUM1=5ETATS<?)*GAMMA*C0MFGA*SUM1 
SUM?=8ETAN*GAMMA*C0MFGA* (-1 « ) »*NS*SUM2 

SUM3=SUM3 ♦ AMACH* <3. 1 41S926 tt NS) **2./2./ ALENGTH*AMU (LS .NX « 1) 

SUM3= (-1 . > »AMACH*SUM3 
DUM1=CMPLX<0. 0.0.0) 

IF (LS.EQ.LHAT) 0UMl=<BETAtS(2) - BETAWI ) *AKl 
OUMl=DUMl +DUM3* <G1 (X2.N5) -GHXI.NS)) 

OUM1 =DUM1/PSI (DUMMY) /F NO°M (LHAT .MHAT .NHAT .RLMOA . ALENGTH) 

AMU (LS .NX .2) = (DUM1 *SUM1 *SUM2 + SUM3 ♦ SUM4) / <0MFGA**2 - ETAd.HAT 

1 .MHA F .NX. ALENGTH .BESPRIM (MHAT* 1 *LS) ) ) /EPS I ZN( LHAT .MHAT *NS. ALENGTH, 
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ZRLMOA) 

BO CONTINUE 

IF < J .EO. ? ) GO TO 14 
IS DO 11 L * 1. LOEX 
00 11 N as 1 • NOEX 

11 AMU (L »N« 1 1 = AMU(L*N»?) 

14 IF ( MATRIX .NE. 3HYFS ) 00 TO 13 
KS = 1 
KF = LOEX 

WRITE (6*100) J* KS* KF 
DO 12 N : 1 , NOFX 

12 WRITE (6. 101 ) N, ( AMU(L.N.l) * L = KS, KF ) 

13 CONTINUE 

IF ( AMU(L0EX,NDEX*1) .Nr. AMU(LDEX,NDEX.2> > GO TO IB 
CALL BETAIJ ( J*1 ) 

IF < ISET .ME. 0 > GO TO 17 
60 CONTINUE 
17 CONTINUE 
JP = J 

100 FORMAT ( *0 N J EQUAL *,I2.» L EQUAL *,I1»* TO *.I2) 

101 FORMAT ( * »,T3.» **10G13.6) 

RETURN 

END 
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FUNCTION E°SIKAT (L»M,N,ALENGTH,RLMOA) 

DUM1 s ALENGTH 

IF ( N .NE. 0) OUMl = OUM1/2. 

IF ( M *EQ. 0) GOTO 1 

EPSIKAT = (0.5-M *M / (?.*RLMOA*PLMDA) ) #m;Ml 
RETURN 

1 EPSIKAT = OUMl/?. 

RETURN 

ENTRY EPSTZN 

IF ( N .EO. 0 1 GO TO 2 

EPSIKAT = ALENGTH / ?. 

RETURN 

? EPSTKAT= ALENGTH 

return 

ENTRY ETA 

EPSIKAT = <N- 1 ) * IN- 1 > * 7.1A15926 *3.1MS9?6 / (^LENGTH * ALFNGTH) 
i*RLMOA * RLMOA 
RETURN 
END 
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COMPLEX FUNCTION GI (X,N) 

REAL LENGTH 

COMPLEX BETA1, BETA?, BETAN, BETAC » BET At, BETAWI. OMEGA, C, A" l . 

I AK? » PSI? 

COMPLEX OUM1, GOIIM, flETA 

COMPLEX Mil , BET4IS ♦ I 

COMPLEX D1 , 02 ♦ 01 , 04 , D5 . 06 

COMMON / RLKA / BFTA1, BFTA2, AMACH. GAMMA , OMEGA. ALENGTH. C, AKl 
1* A < 2 * BETAN, PLMDA, BET A I , BETAC, BETAWI, LHAT* MH4T • NHAT, 

1 PSI2 » XI, X? 

COMMON / BLKB / MU(1,30.p> , BETA IB ( 2> , I « MATRIX 
COMMON / RLKC / LENGTH , 01 (TO) ♦ 0?(1C> , 01(30) f 04(10) . OB(30 
$) » 06(10) ♦ Al(?) 

GDUM(A*I,BETA.X)=(CEXP(I*RETA*CMPLX(X,0.O) ) * ( I*RETA*CMRLX (COS (A*X > 
$,0.0) ♦CMPLX(A«SIN(A*X) ,0.0) ) ) / (CMPLX (A*A.O.O) -BET A* BET A ) 

A = N * 3. I4lS9?653Sn979 / ALENGTH 

Gl = <0MEGA»AMACH*BETA1)»GDUM<A,I,8ETA1,X) 

IF ( N ,E0. 0 .ANO. CARS < BETA2) .LE. 1. E-100 ) GO TO 1 
Gl = I*GAMMAMGIM0MEGA,AMACH*HETA2)*C*GDUM{A, I.RETA?,X> ) 

return 

1 BETA = C * OMEGA * X 

Gl = GAMMA * ( Gl ♦ BETA ) * I 
RETURN 

ENTPY G1FBST 

A = N * 3. 141S926B3SB979 / LENGTH 

IF ( N .EO. 0 .ANO. CABS ( BETA2) .LE, 1. E-100 > GO TO 5 
Gl = < CAMACH*4MACH-I. > *RFTA 1 *9ETA K?. *OMFGA*AM ACH*BETA 1 +OMEGA»OMFG 
5A -MHAT*MHAT)*GDUM(A,I,8ETA1+BETA1 ,X) ♦ ( ( AMACH*AMACH «l.)* 

$2.*C*8ETA1*8ETA2*2.*0MEGA*AMACH»C*(RETAKBETA2) ♦ (OMEGA*OMEGA-MHa T® 
$MHAT)«?.*C)*GDUM(A,I,RETA1+BETA2.X> ♦ ( ( AM ACH*AM ACM -1. 

$)*C*C*BETA2*RETA?*2.*0M£GA*AMACH*8ETA?*C*C+ (0MEGA»0MEGA-MHAT*MHAT) 
$*C*C> «GOUM (A, I ,flETA?*BET42,X) 

RETURN 

<5 Gl = ( (AMACH*AMACH-l.>*BETAl*8ETAK?.#0MEGA*AMACH«BETAK0MEGA«OMEG 
$A -MHaT*MH4T)*G0UM(A,I.RETAKBETA1,X) ,{( AMACH* AMACH -I.)* 

$2 »*C*flET A 1 *8ET42+2, *OM£GA* AMACH*C* ( BETA 1 ♦BET A2) ♦ (OMEGA*OMFGA -MH 

SMHAT.) *2.*C»-»GDUM (A, I , BETA 1 *BET A2 , X ) ♦ (OM£GA*OMEGA“MH4T*MHA 

$T)*C*C*X 
RETURN 
ENO 
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COMPLEX FUNCTION EPSTRZN (L*M * N» ALENGTH ♦ RLMOA > 
COMPLEX OUM? 

OUMl = i. 

if(m.eq.O) r>UMi = 

0UM?=FNORM<L»m,N»9LmDA. ALENGTH) ** 2 
EPSlR/M=D(jMl * OUM? / 3.14159265358979 
RETURN 
END 


complex function fnopm (l n . rlmd a » alength) 

L)UM 1 = 3*141 59?6 

IF ( M .EG. 0 > OUMl = 6.P831852 

Ol.!Mi=DUMl*EPSIKAT <L *M ,N ♦ aLENGTH ♦ RLMna ) *6E5SEL (M+l «L) **2 
FNORM = CSQRT 1 CMPLX ( OUMl * 0.0 ) ) 

RETURN 

END 


FUNCTION F0(X) 

F0=. 79780454-. 000 00077* < 3. /X>-. 00552740* (3. /X>**2 

1- .000095I2* (3./X>**3*.001 37237* <3./X>**4 

6 4072805* ( 3. /X ) »*5*. 00014476* (3. /X>**6 
RETURN 
ENT D Y FI 

F0=. 79788454*. 0QOO0 154* < 3. /X) *.01659447* <3. /X)**? 

1 + .0 00I710S* (3./X) **3-.00P495U* (3 . / X ) **4 ♦ . 00 1 1 3653* ( 3 . /X 1 **5 

2- .00020033*(3-/X»**4 
RETURN 

END 
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COMPLEX FUNCTION G2(N> 

REAL LENGTH 

COMPLEX BETAI* BETA?* BET AN » BETAC* BETAI » BETAWt ♦ OMEGA* C* AK 1 » 
l AK?* PSI? 

COMPLEX MU * AMP * I 

COMPLEX D1 ♦ 0? ♦ D3 ♦ DA * 05 * D«S 

COMPLEX GRUM * BETA 

COMMON / RLKA / RETA1* 8PTA2* AMACH* GAMMA* OMEGA* ALENGTH* C« AK 1 
1* AK2* BETAN* 9L«0A* BET A I » BETAC, RETAW I « LHAT* M HAT • NHA T » 
i PST2*Xi* X? 

COMMON / BLKR / MU(3*30*?> * AMP( 2) *1 * MATRIX 

COMMON / aLKC / LENGTH « 01(30) ♦ 02(30) ♦ 03 (30) * . 04 < 30) « OR(30 
$) ♦ D6(30) * A 1 ( ? ) 

G0UM(A.RETA.ALENGTH*N> = (BETA*(CEXP(BETA*CMPLX (ALENGTH* 0. 0> )*CMPt_X ( 
$-1 *0*0*0) (N*?) +CMPLX (-1 *0*0« 0> > ) / (CMPLX (A* 0*0) *BETA*BETA) 

A = N*N*3. 1A1S92653SRR7P*3. 1415926S3BBP79/ALENGTH/ALENGTH 

G? = -BETAI* ( AMACH«BFTA1 ♦?.*OMEGA) *G8UM ( A.BETAl * I , ALENGTH *N > 

IF ( N *E0* 0 -ANO. C ABS ( 8ETA2) *LE. l.E-10 ) RETURN 

G2 s G2-C*dETA?*(AMACH*BETA2>2.*0MEGA) *G8UM(A*BETA2*I*ALENGTH*N) 

RETURN 

ENTRY PSI 

A = NHAT*NHAT*3.lAlSQ?t>53SB97P*3. I A1SQ2B5358979/ ALENGTH/ ALENGTH 

G? = GBUM ( A . I * BETAI * ALENGTH , NHAT ) 

IF ( NHAT *E r ). 0 .ANO. CARS( RETA2) .LE. 1. E-100) GO TO 3 
G? = (G2*C*GBUM (A*I*RETA?» ALENGTH *NHAT) )/EPSIZN(LHAT .MHAT .NHAT .ALE 
$NGTH*RLMDA) 

RETURN 

3 g? = ( G? + C * ALENGTH ) / EPSI ZN ( LHAT . MHAT «NH A T . ALENGTH . PLMO A > 
RETURN 
END 
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REAL FUNCTION INTEGRL <F,A*8,LHAT»L*L8AR,MHAT,M9AR> 

oimfnsion W ( 1 | > * X ( 1 1 ) 

OAT A (W ( I) ,1=1, 11 >/.02l?9lO 1837554,, 05761665831 124«.09340039B?7R2 
15. . 12052016961432, . 13642490095628. . 141 49370«92875» . 1364?4900Q56?a . 

2.12352016961432. . 09340039827825.. 05761665831 1 24.. 021 291 01 837554/ 
DATA <X( I) ,1 = 1, 11 >/.007o973l9952S9.. 04691007703067,. 1 2291663671 4S 

1 9.. ?3C765344O47l6,.3601847934l911..900P0000000000.. 63981 520659089. 

2.76923465505284.. 87708336328542. .95308992296933.. 99204268004742/ 
CONST = 8 - A 

ARGl = 6ESPRIN (MHAT ,L H AT) 

ARG? = BESPRTM(MHAT,L> 

ARG3 = 8£SPRIM(M8AR,L8AR) 

INTFGRL « W( 1) *F(A*CONST*X( 1 > , ARC, 1 , 4RG2 , ARG3 > 

00 1 I = 2,11 

1 INTFGRL = INTEGRL ♦ W( T)*F(A+CONST*X ( I > , ARGl ,ARG2,AR63> 

RFTilRN 

EMO 
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HEAL FUNCTION 11 ( X . RLMDA ♦ RLMOAL * RLMRAR ) 

C *** 

c *** THIS FUNCTION SURPROGPAM ASSUMES THAT MRAR EQUALS TWO MHAT 
C *** AMO that MHAT EQUALS ONE 
c *** 

REAL JO * Jl * JR 

J2( ARG ) = ?.*J1 (ARG)/ARG-JO<ARG) 

C 

c *** these are THE integrals for m equal zero 
c *** 

Il=Jl(RLMOA*X)*JH.RLMDAL*X> tt JO(PLMBAR<»X>/X 4 RETURN 
ENTRY 12 * IlsJl (RLMOA»X>*Jl (RLMOAL*X)«-JO(PLM8AR*X)*X 4 RETURN 
ENTRY 13 4 U = JO(RLMr»AL*X)*Jl (RLMOA»X>*JO(RLMBAR*X) 4 RETURN 
ENTRY 14 4 I1=JO(RLMOA*X)*J1 (RLMOAL*X>»JO<RLMBAR»X) 4 RETURN 
ENTRY IS 4 I1=X*JO(RLMOA*X)*JO (RLMDAL»X)*J0(PLMBAP*X) 4 RETURN 

C *»■» 

c these are THF INTEGRALS FOR m equal two 
c *** 

ENTRY lb 4 I l=Jl (RLMOA*X) *J1 (RLMOAL*X) W J2 (PLMBAR*X) /X 4 RETURN 
ENTRY 17 4 I1=J1 (RLMOA*X)*Jl (RLMOAL»X)*J2(RLMBAR»X)»X 4 RETURN 
ENTRY IB 4 Il=J0 (RLMnAL*x )*J1 (RLM0A*X)*J2(RLMBAR*X) * RETURN 
ENTRY I q 4 H=JO(RLMnA*X) *J1 <RLMDAL*X)»J?(RLMBAR»X) 4 RETURN 
ENTRY HO 4 Il=X*JO (RLMOA*X) *J0 <RLMDAL*X) *J? (RLMRAR*X) $ return 
ENTRY 111 4 II = X*Jf>(RLMr)A*X)*JO{RLMBAR“”' >J1 (RLMOA*X) 4 RETURN 

ENTRY II? 4 II = JO (RLMR*R*X> *J1 (PLMQA»X) **2 4 RETURN 

ENTRY 113 S II = X*JO(RLMQA*X)*Jl (RLN0A*X)*J2<RLNBAR*X) 4 RETURN 

ENTRY 114 4 II = Jl <RLMDA*X)**2*J2<RLMBAR*X> 4 RETURN 

ENTRY 115 4 II = Jl (RLMDA*X>*»2*J1 (RLMBAR»X> 4 RETURN 

ENTRY Rl 4 II s X* ( JO (RL M OA*X ) * Jl ( RLMQA*X ) ) **2 4 RETURN 

ENTRY R? 411= JO (RLMf)A*X) *J1 (RLMDA*X> **3 4 RETURN 

ENTRY R3 4 II = Jl <PLMOA*X > **4/X 4 RETURN 

ENTRY R4 4 II * Jl (RLMOA*X > »*4*X S RETURN 

ENJOY RS 4 II = JO(PLHOA«X)*«3*J1 (RLMDA*X) 4 RETURN 

ENTRY R6 4 II » ( JO ( RLMDA*X > * Jl (RLMr)A*X ) > **2/X 4 RETURN 

ENTRY P7 4 I l = JOJRLMOA*X) * Jl <RLMOA*X) **3/X/X $ RETURN 

ENTRY RS 4 11 s Jl <RLMt>A«>X) #*4/X/X/X $ RETURN 

ENTRY R9 4 II = ( 1 . -PLMR AR*RLMB AR/RLMDA/RLMOA ) *BFSSEL ( IF I X < RL'-tq AP* 
4 1 . ) <TEI X (RLMOAL) ) **?/2, $ RETURN 

ENO 
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REAL FUNCTION J0<X> 

IF (X-3. )50*S0*5l 

*50 JO*l--2.269QQQ7*CX/l.)**?*l*2656?08*CX/3.)**4-.316l8<.fi*(X/3.»**ft 
l*.044447O*(X/3.)**8-..003Q44A*(X/3.) ** 1 0* . 00021 0 0* (X/3. ) **1? 

RFTURN 

51 J0=(F0 (X)*CO5(TH0(X) > >/SOPT(X> 

RETURN 

ENT3Y J1 

IF(X-3.)S?»S?.53 

52 J0= <.5-.Sft24RQ85* (X/3. )»* 2 *. 21 093573* ( X/3. ) **4- . 039S428R* < X/3.) 

• 004433 1R* ( X/3. 00031 761* ( X/3.) **10*. 0000 1109«( X/3. ) ** l? 
2 ) *X v 

return 

53 J0=(F1 (X)*COS(THl ( X ) > ) /SORT ( X ) 

return 

END 
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/ 

SUBROUTINE LINEUP ( LENGN’O * F * G ♦ AK * FIMC . FMAX * JX ) 
COM a LEX BETA 1 « BF.TA?, BE T AN ♦ BET AC ♦ BETAI* RETAWl* OMEGA. C* AK1. 

1 AK?* PSI?. AMU* BETAIS. C12 

COMPLEX G* AN* AK* J* AA * A5* BETA* OtJMl* 0UM2, PSI* OUMMY. Gl. 

1 02 *F * AKA .FINC* FMAX 

COMMON / QLK A / BETA!* RFTA?* AMACH. GAMMA* OMEGA* ALENGTH* C* AK 1 
1* AK2* HETAN* PLMOA* BETAI* 8ETAC* RFTAWI* LHAT* MHaT* NhaT. 
1PSI?*XI*X? 

COMMON / 8LKB / AMU ( 3*30*2) t BET Al S ( 2)»I»MATPIX 
COMMON / BLOCK / C12<30*30) * ISET*L0FX*N0FX*FRQ0Ri.FWR0P2 
I = CMPLX ( 0.0 » 1.0 ) 

RLMOA = BESPPIM ( MHAT ♦ 1 * LHAT ) 

NHAI = NHAT ♦ l 
OMEGA=F*l. BA 11*378 

BETA=CSORT ( OMEG A*AMACH*OMEGA* AMACH* ( AM ACH*AMACH- ] . ) * (SLMDA**? - 
10MEGA*0M£GA)> /{\ .-AMACH* AM ACH) 

BETA 1= (OMEGA*AMACH>/ (1.-aMACH*AMACH) .BETA 
BET A2= (BETA 1 > -2. *RETA 

C * - CEXP ( I * ALENGTH* ( 3ETA1 - BETA2 ) ) 

C = C * (BETA! ♦ BETAN * GAMMA * (OMEGA ♦ AMACH * BETA 1 >) / ( 

1 BETA? ♦ BETAN * GAMMA * (OMEGA ♦ AMACH * BETA2 ) ) 

PSlP = PSH DUMMY ) **2 

AKi=(OMEGA* AMACH *RETA1+C*(0MEGA* AMACH*BETA?))*GAMMa 

AK 1 = CMPLX(-AIMAG(AK1) *RFAL(AK1) ) 

DUMl = CMPLX(-AIMAG(BETAl ) .REAL(BETAl)) 

DUM1 = 0(JM1*ALENGTH 

DUM? s CMPLX (-AIM AG (RE TAP) » REAL (BETA 2 ) ) 

DUMP = OUM 2* ALENGTH 

AK2 st ( (0MEGA*AMACH*BETA1 >*CEXP(DUM1) ♦ (OMEGA+AMACH*BETA2> *C*CEx a ( 
IOUM?) )*GAMMA 

AK2 = CMPLX(-AIMAG(AK?) ,REAL(AK2) ) 

A4 = AK1/EPSIZN(LHAT.MHAT*NHAT*ALENGTH,RLM0A)/PST (OUMMY) 

ARsRETAC* (G l (X2. NHAT) -Gl (XI *NHAT) ) /EPS IK AT (LHAT, MHAT.NHAT* ALENGTH* 
lRLMOA) /PSI (DUMMY) 

BETAWI = (OMEGA**? - ETA (LHAT, MHAT.NHAI « ALENGTH. PLMOA) - (AMACH*G2 ( 

1 NH AT > /EPSIZN (LHAT .MHAT.NH AT .ALENGTH, RLMOA) + RET AN* AKP/EPS I ZN (LHAT * 
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1 MH AT ♦NHAT»ALENGTH*RLMDA) * (-1 . > ** (NHAT*2> > /PS I I DUMMY) > /A4 
BETAl = BETAWI - AS / A4 
BET A IS ( 2 ) s BETAIS ( 1 ) = 3ETAI 
WRITE (6, 100) LENGNO 
WRI TE (6 • II 0) LDEX ♦ NOEX 
WRITE <6* 1 16) AMACH 
WRITE (6« 117) GAMMA 
WRITE (6. 118) OMEGA 
WRITE (8. 11R) ALENGTH 
WR I TE (6*101) LHAT * MHA T « NHAT 
WRITE <6*1 15) X 1 * X2 
WRITE (6* 1 13) F 
WRITE (6« 102) 8ETAN 
WRI TE (6 • 1 20) G 
WRITE (6*121) 8ETAC 
WRI TE (6 • 108 > AK 
WRITE (6» 105) BETA W I 
AN * 1. / GAMMA - BETAWI / AMACH 
WRI TE (6 ♦ 1 03) AN 

6FT4WI = ( BETA1 * C. « BFTA2 ) / ( GAMMA * (OMEGA * ( l. ♦ C ) * 

1 AMACH * ( 9ETA1 ♦ C * BFTA2 ) ) ) 

WRITE(6*U2) BETAWI 

AN a 1. / GAMMA - BETAWI / AMACH 

WR I TE (6 * 1 03) AN 

WRI TE (6 * 109) BETAIS(l) 

AN = 1 • / GAMMA - BETA IS ( 1 ) / AMACH 
WRITE (6* 103) AN 

IF ( X?.LE.1.F-10.0R.CABS(BETAC) .LE.l.E-10) GO TO 1 
JY = JX - 1 
CALL CALMU( JY ) 

1 JX = JY «■ 1 

IF ( JX ,EO* 30 ) 2 * 3 

2 WRI TE (6* 1 22 ) 

3 RETURN 

100 FORMAT(*I THIS IS ENGINE NUMBER *A1 0) 

101 FORMAT (*0 THF PRIMARY MODE ASSUMED IS LHAT = **!?«* MHAT = **!?»* 



D 





INHAT = *.I2) 

102 FORMAT (*0 8FTAN = **2G21.14> 

103 FORMAT <*♦*» 70X.* M = * * 2021 . 1 4. * I*> 

105 FORmAT(*i) 8FTAWI = *«?G?1.14> 

107 FORMAT (*0*.G2J . 14.* WAS THE TIME FOR THE ITERATION *.021.14) 

108 FORMAT (*+*.70X»* K = *.2021.14) 

109 FORMAT <*0 8FT A I ( 1) = *.2G?1.|4> 

110 FORMAT (*+*,70X.* THE COEFFICIENT MATRIX IS *«I2.* X*,I3) 

112 FORmAT(* 0 BETAWIO = *.2G2l.l4) 

113 FORMAT (*0 W/WO = *.2021.14) 

115 FORMAT <*+*.7PX«* THE LINFR BEGINS AT *.F6,4.* ANO ENOS AT *.F6.4) 

116 FORMAT (*0 THE MACH NUMBER IS *.G21.14> 

117 FORMAT (*+*.70X«* THE RATIO OF SPECIFIC HEATS IS *.G2l.l4) 
lift FORMAT («0 THE FREQUENCY TS *. 2G21.14) 

]19 FORMAT (*+*»70X«* THE LENGTH OF THE COMBUSTOR IS *.G21,14) 

120 FORMAT (* + **7!)X.* G = *. 2G21.14) 

121 FORMAT (*0 BETAC = *, 2G21.14) 

1 22 FORMAT ( *1 * 1 35 ( 1H* ) /* *47(1H*)* THE COMBUSTION RESPONSE DIO NOT CON 
SVERGE *46 ( 1 H* ) /* *l3^<lH*n 

END 


FUNCTION THO ( X > . 

TH0=X-» 785398 1 6-. 04 166397* ( 3 . /X )-. 09003954* C3./X) **?♦. 00262573* 
1 (3./X) **3-. 00054125* (3./X) **4- . 00029333* < 3. /X ) **5+ . 000 1 3558* 

2 (3./X) **6 
RFTURN 
ENTRY TH1 

TH0=X-?. 3S6 19449*. 1249961 2* (3. /X) +.000 05650* (3. /X>**2 

1- . 00637879* (3. /X)**3 + . 0007434 8*(3./X) **4 + . 000798 24* < 3. /X > **5 

2- .0A029l66*(3./X)**6 
RETURN 

END 
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SUBROUTINE CONST 

REAL II, I?, 13. IA« IS, 16.17, 18. 19, II O.INTEGRL. LENGTH 
COMPLEX Dl , 02 . 03 , D4 , 05 ♦ 06 
COMPLEX A , R . 0 » 1 516 < 10 .3 ♦ 2) ♦ GlFRST 

COMPLEX HI. H?, I • 6ETAN. BETAC. RETAI. 8ETAW I * OMEGA. C. AK { « 

1 AK?» PSIR. MU . AMP , CI2.EN0RM 

COMMON / RLKa / R1 • RP . AMACH. GAMMA. OMEGA. ALENGTH, C« AK l 
l* A<2» HETAN. RLMOA , BETaI. BETAC. RETAWI. LHftT, MHaT. NhaT. 
IPSIR.XI.X? 

COMMON / RLKfl / MIJ(3.30.P> ♦ AMP( P) .1 . MATRIX 

COMMON / PLKC / LENGTH , 01(30) . 02(30) . 03(30) . 04(,in) , 0s <30 
*) , 06(30) . A 1 (?) 

COMMON / BLOCK / Cl? (30 .TO) *ISET « LOE X . NOE X. ERROR 1 .ERROR? 
EQUIVALENCE ( 01 ( 1 ) , 1516 ( 1 .1 , 1 ) ) 

EXTERNAL II . IP.I3.IA.I5. T6.I7, l8.T9.riO 

D(A.B.X.N) = (A*(CEx p (A*CMPLX (X.O.O))*CMPLX (-1 . 0 . 0 . 0 ) ** (N+P) .CMPLX 
I (-1.0. 0.0) > )/(A*A.8*R) 

C *** 

C **» NOTE THAT IN THIS SUBROUTINE BETAI . BETA? BECOME 81 . B2 
C *** RESPECTIVELY. 

C *** 

03 * HI ♦ B? 

00 1 NX = 1 , NOEX 
N = NX - 1 

B = CMPLX ( N * 3.141592*5358979 / LENGTH . 0.0 ) 

01 (NX)=B1*0(?.*I*8I,R,LENGTH,N) *B3*C *D (I *R3.B .LENGTH, N) *C «C 
$ *R2*D(I*2.«B2,B, LENGTH, N) 

D? (NX> =81*R1*R1*0(2.*I*R1 .8, LENGTH. N) .B i*R2*B3*C *0 ( I *B3 .B. LFNGT 

*H,N> ♦B?«8P*BP«C *C *r>(2.*I*R2, Q, LENGTH. N) 

D3(NX)=0(?. 8. LFNGTH.N) ♦C«?.*0( I*B3.B, LENGTH. M) ,C *C *D 1 1 
$*?.*«?, B.LENGTH.N) 

OA (NX)=B1*81*D(?.»H>R1 . R, LENGTH, N) *C* (81 *BI +R?*B?» *0 < 1*R3.B .LENGTH 
S.N).C»C *B2*R2*D<2. *1*82, R, LENGTH, N> 
l OS (NX)=81#bi«0(?.*I*Rl .B, LENGTH, N) *2,*C *R1*B?»0( I*R3 ,R, LENGTH ,N 

*>♦0 *C *B2*B?*0(P. *1*82, R, LENGTH, N) 

IF ( CABS ( B? ) .LE. l.E-100 ) A , 5 
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4 Ol( 1)*81*0<2.*1*BI,8,LENGTH.O) *B3*C *0 ( 1*83. 8. LENGTH. 0 ) 

D2( l)sBl*Bl*Rl«0(2.*I*8! ♦ B, LENGTH. 0) 

03 ( 1 > =0<2.*I*B1,B, LENGTH, 0 > *C*2, *0 ( I*B3 . B , LENGTH. 0 ) ♦C*C*LENGTw 
04 ( 1)=R1*B1*0(?.*I*B1 »B, LENGTH .0) ♦ C*(B1*BKB2*R?>*0 ( I *83*3. LENGTH 
S .0) 

08 ( 1) =8 1*8 1*0 (2. *1*81 *8. LENGTH, 0) 

5 00 2 L = 1 . LOEX 

MU (L * 1 1 ♦ 2 ) = CMPLX ( I NTEGPL ( 1 1 * 0 » 0 ♦ 1 • 0 » LHAT *LHAT ,L *2, l ) . 0.0> 

MU (L* 12*21 = CMPLXdNTEGOL (12.0. 0.1. O.LHflT.LHAT.L.2.1) » 0.0) 

MU (L ♦ 1 3 * 2) = CMOLXUNTEGDL (13*0. 0.1.0. LHAT.LHAT.L. 2.1) . 0.0) 

MU(L * 15*2) = CMPLX ( I NTEGPL (15*0*0. 1.0.LHAT»LHAT»L»?.l) * 0.0) 
MU(L*16,2> = CMPLX(INTEGOL(IG.0.0*1.0.LHAT*LHAT.L.2.3) . 0.0) 
MU(L*17.2) = CMPLX ( I NTEGPL (17*0.0*1.0 . LH AT, L HA T»l«2, 3) . 0.0) 
MU(L*1B.2) = CMPLX<INTEGPL(I8.0.0.1.0.LHAT,LHAT.L.2*3) . 0.0) 

2 MU(L.20.2> = CMPLX(INTEGPL(I10*0.0.1 .O.LHAT, LHAT.L, 2.3) . 0.0) 

: c 

; C *** CALCULATION OF INTEGRALS TS AND 16 . ANO THEN StJM 
C 

I NBA I = NRAR ♦ I 

MU (2*1,2) = CEXP ( B 1* I *LENGTH) 

MU (3*3.2) » CEXP(B2*I*LEMGTH) 

MUU,t,2)=CMPLX(l.O*0.)/PSI2/FNORM(LHAT,MHAT,NHAT,RLMOA.LENGTH)o*2 
MU (2.2,2) * BETA I* ( 1 * *C) **2 

MU ( 1*4.2) = 8ETAN*(WU(2,1,2MC*MU(3.3,2> )**2 

MU (2*3.2) = BETA I* ( <RKC*R2) **2* ( AM ACH* AM ACH-1 . ) +2. *OMEGA* AM ACH* { 1 
S. *C) * (B1 *0*82) +OMEGA*OMFGA* ( 1 . *C) * (1 . *C> ) 

MU( 2*4*2)= BET AN* ( (R1 *MO (2, 1.2) +82*C*MU ( 3, 3 . 2> ) **2* ( AMACH*AMACH- 1 . 
$) ♦2.*0MEGA*AMACH*(MU<2,1 .2) ♦C*MU(3*3,2) > * (B1*MU (2,1 .2) .B2*C*MU(3.3 
*♦?> ) ♦0MEGA*OMFGA*(MU(?.l .2) .C»MU(3.3,2) ) **2) 

MU (3*1*2) = BESSEL (MHAT+ 1 .LHAT) **2 * BET AC 
DO 3 NX - I , NDEX 
N = NX - J 

MlJ ( l *2*2) = ( l .—GAMMA ) /4 * * (AMACH*AMACH*AMACH*02(NX) +2.*OMEGA*AMA 
$CH*AMACH*(DS(NX) *D4(NX) /?.) ♦3.*0MEGA*OM£GA*AMACH*01 (NX ) +0MEGA**3 
$ *03 (NX ) ) - (AMACH*D2(NX)*OMEGA*05(NX) )/2. 

MU < 1 » 3*2) = (AMACH*Ol (NX ) * OMEGA * D3(NX>>/2. 
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MU(3 9 2*2> - GIFRST't X? « N 5 - G1FRST ( KI ♦ M I 

>iun»*»?) = cmplx ( ( -i. ) ** n ♦ o.n > 

MU < 1 *5 1 2 ) = MU(2»2«2)*MUng4«2>*MU(3,A.2) 

MU(?*S*2> * MU (2 « 3»2> + MU f ?♦ 4» 2> *MU ( 3*4» ?> 

OO 1L s 1 I LOEX 

I 51 A (NX *L* 1) =A1 ( 1 >*MU(1 «l »2 ) /FNOPM (L * 0 * N tBESPRlM ( 1 ,L) *LEMGTH) * ( 
SMUIL* 12*2) .2#?) ♦3»2>*( p LMDfl*|OLMDA*MUCL.lS.?> 

j- 2 .*MHAT*MU(L» 1 3# 2) ) ♦CMPf.X(MHAT*MHAT*2. • 0. 0 ) *MU (L • 1 1 ♦ 2) > 
S+GAMMA/S** (MU( l.S*2>* (RLMOA* (RLMDA*MU(L*15«2) -2 . «MhAT«MU (L ♦ 13 »?) > 
{♦CMPLX (MHAT*MM4T*2. ,0.0) *M(J (Lt 1 1*2) ) *MU (L*12.2)*MU<2«5.?> 
?*MU<3»1»2) *B£SSEL ( l *L) *MU (3*2*2) ) ) 

3 I5IMNX*L.2)*A1 (2) * MU(}*1*2> / FNOPM (L . 2 *N .BESPPIM ( 3 . L > tLENOTH) 
$* (MU(L*I7*2>*I*MM(1*?*2> ♦I*MU(1 .3*2) «PLMDA* (»LMDA*MU(L« 20.2) 

$-?.*MHAT*MU<L.lP*?> ) ♦GAMMA/P,* (MU (1 .5.2) *PLMOA* (HLM0AO 

$MU(L*?0*2>-?.*MHAT*MU(L*tft»2>) ♦MU(L*17*2)*MU(2*5»?> 

?*MU (3*1*2) *8ESSEL (3»L> *MU <3* 2*2) ) ) 

RETURN 

END 
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COMPLEX FUNCTION GAMM < L ♦ M , N > 

REAL LENGTH 
COMPLEX F OF X . I 

COMPLEX BETA 1 * BETA?, BETAN, BETAC, BET A I * BET AWT * OMEGA, C, AK l , 

1 AK? » PSI?» MU , AMP , Cl? 

COMMON / BLKA / BETA 1 « B C ’TA2* AMACH, GAMMA. OMEGA* LENGTH, C* AKI 
1, A<2, BETAN, RLMQA » BETaT, 9ETAC, PETAWI * LhaT, MHAT, NHAT. 
1PSI?,XI ,X? 

COMMON / BLKR / MU < 3 . 30 * 2 > * AMP ( . ? > , I , MATRIX 

COMMON / BLOCK / Cl?(30*10) *I5ET,LDEX*N0FX,FPR0P1*ERR0R? 

GAMM - AMACH*AMACH*N*N*3. 141592653SA979*3. 141592*5358979*3. 1415926 
?G35«979/LENGTH/FN0RM(L,M.N,8ESPRIM(M*1 ,L) .LENGTH) **?*RESSEL ( M ♦ 1 ,L ) 
5**2 

IF ( M > ? , i , 2 

1 RETURN 

2 GAMM = GAMM* (1 «-M*M/BESPPlM (M*l ,L) »*2) /2. $ RETURN 
ENTRY DELTA 

GAMM =_ -GAMMA*I*?.*0mEGA* 3. 1 4159265358979*8ESSEL (M* 1 ,L ) **2/FN0RM ( L 
$,M*M*8ESPRI M fM*I,L) , LENGTH) »*? 

IF < N ) 4 , 3 * A 

3 F OF X j CMPLX ( X2 - XI , 0.0 ) 5 GO TO 5 

4 F OF X = CMPLX (0.5,0. 0 >*< X2-X l ♦LENGTH/? . /N/3. 14 1 59265358979* ( SIN < N 
$*?.*3.14159?6S358979**2/1 ENGTH) -SIN<2.*N*3. 14159265358979*X1/LFNGT 
SH) ) -AM ACH/I/?. /OMEGA* (SIN (N*3 . 1 4 1 S9265358979*X2/LENGTH ) **2-SI N (N*3 
t. 141 5926535«979*X1 /LENGTH )»*2) ) 

5 IF ( M | 7 , <S , 7 

6 GAMM = GAMM«(BETAI-»BETAN*(-1.)**N+BETAC*2.*F OF X) $ RETURN 

7 GAMM = GAMM*< (BETAI .BETAN* (-1 •> **N) * (1 .-M*M/BESPRI M (M ♦ 1 ,L )**?)/?. ♦ 
SdETAC*F OF X) S RETURN 

ENO 
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SUBROUTINE INIT (A*H«C«L*N) 
DIMENSION A ( 1 ) , 8(1) , CU) 

INOFX = L*N*? 

DO l I = l • INDEX 

1 B < I ) = 1777 0000 0000 0000 OOOOB 
INDEX = INDEX * 2 

00 ? I = I * INDEX 

2 C(I> = 1777 OftOO OOOO 0000 00008 
INDEX = N»N*2 

00 3 I = 1 • INDEX 

3 A < I ) = 0.0 
RETURN 
END 
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subroutine nonline ( locate > 

REAL LENGTH 

COMPLEX GAMM , DELTA * MUSNO. MUCROSS ( 3* 30 . 2> , FNQRM* COR 
COMPLEX 8ETA1 ♦ BETA?* BETAN* 3ETAC » BETA I * BETAWT • OMEGA. C. A*l. 
l AK?t PSI?* AMU* BET A IS * C12* I* ETA 
COMMON / fiLKA / BET A l * BETA?* AMACH, GAMMA. OMEGA* ALENGTH, C» AK 1 
1* AK2* BETAN, PLMDA, BETaI* BETAC* BETAWl. LHAT* MHAT. NHAT. 

ipsi?*xi*x? 

COMMON / BLKB / AMU ( 3*30*2) .BETfllS ( 2) *1, MATRIX 
COMMON / RL«C / LENGTH , MUSND (30 *3*2) * A 1 C 2) 

COMMON / BLOCK / CI2(30. 30) *ISET*LDEX,NDEX.ERR0R1 .ERROR? 
EQUIVALENCE (Cl ? ( 1 * 1 )* MUCROSS < 1 * 1 * 1 ) ) 

C *** 

C *** CALCULATE NECESSARY CONSTANTS 
C *** 

IE ( LOCATE oEQ* 10HCROSS PROD ) CALL X PROO 
DO 1 NX = 1 o NDEX 
N = NX - 1 
00 1 L = 1 , LOEX 

COR = A.*OMEGA*OMEGA-ETA (L* 0 ♦ NX* LENGTH* BESPRIM (1 *L) ) 

MUSND (NX * L * S > =( MuSND ( NX *L ♦ 1 > ♦ MUCROSS (L *NX . 1 ) ) / ( COR + GAMM (L«0 
S.N) ♦ DELTA(LoO*N) ) / EMORM (L *0 *N. BESPRIM ( 1 ,L) .LENGTH) 

COR = A.»0MEGA#0MEGA-ETA(L*?*NX*LENGTH,BESPRIM(3,L>) 

1 MUSND (NX *L ♦ ?) =(MUSNO(NX,L*2) ♦ MUCROSS (L .NX . 2>> / ( COR ♦ GAMM (L » ? 
?*N> * DELT A (L o ? « N> ) / FNORM ( L * 2 »N * BESPR IM ( 3* L ) * LENGTH) 

IF (MATR I X*NE o 3HYES) GO TO A 
WR I TE ( 6 *60' : ) LOEX 
DO ? NX = 1 . NOEX 
N = NX - l 

2 WRITE (6*601 ) N * (MUSND (NX *L* 1 ) *L=1 »LOEX) 

WR I TE ( 6*602 ) LOEX 

00 3 NX = 1 , NDEX 
N = NX - 1 

3 WRI TE (6 * 60 1 > N 9 (MUSNO(NX,L*2) *L=1*LDEX) 

4 RETURN 

600 FORMAT <*0 N*1?X *L EQUAL 1 TO *I2*6X* M EQUAL 0 ») 

601 FORMAT(* *13* *1 0G13.6) 

602 FORMAT (*P N*1?X *L EQUAL 1 TO *12. 6X* M EQUAL 2 *> 

END 
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subroutine x proo 

REAL I1.I?*I3.14*I5*I6»I7.IB.I9.I10« I NTEGRL « LENGTH 
COMPLEX MU''ROSS(3.30.2> 

COMPLEX MUSNO ♦ ZCR0SS1 . ZCR0SS2 ♦ PSI 

COMPLEX B 1 * B? » I* 8ETAN * BET AC » RETAI. RETAWI, OMEGA. C. A*l. 

I AK2. PSI ? » MU ♦ AMP . C I ? t FNORM 

COMMON / BLKA / 81 * R2 • AMACH. GAMMA * OMEGA. ALENGTH, c. AK1 

I. AK2. BETAN, PLMDA « SET A I * BETAC. RETAWI • LhaT. mhAT. NhaT. 
IPSI2.X1.X2 

COMMON / PLKR / MU(3.30.P> t AMP( ?) . I , MATRIX 
COMMON / RLKC / LENGTH * MUSNO (30 . 3.2) . A 1 (2> 

COMMON / BLOCK / C12H0.R0) .ISET.LOEX.NOFX.EPRORl .ERROR? 
EQUIVALENCE (02(1*1) ♦ MUCROSS (1.1.1) ) 

EXTERNAL I i. 12. 13, U. 15. 16.17.18.19.110 

C %ss 

C sss THESE ARE THE CONSTANTS TO BE USED IN THIS PART OF THE PROGRAM 
C $*S 

WP I T£ (6*600 > 

MU ( 1 ♦ 1 *2) =CMPLX ( 1 .0.0.0) /PSI (0) /FNORM (LHAT* MHAT .NHAT. RLMO A .LENGTH) 
Ml) (2*1.2) = CMPLX( (1.0-GAMMA) /4.0 . 0.0 ) 

MU (O.l.p) * I * OMEGA 

Mu(l. 2.?) = AMACH*AMACH*R1*81*2.*OMFGA*AMACH#81.0M£GA*OMEGA 
mu (2*2.2) = (AMACH*AMACH*B2*B2*2**0MEGA*AMACH*B?.0MFGA*0MEGA) *c 
MU (3.4*2) = CMPLX ( 3. 141R9265358979 / LENGTH . 0.0 ) 

MU( 3.2.2) = AMACH*MU(3.4,2) 

MU( 1 * 15*2) = MUn. 2 . 2 ) * MU(3.2.2) 

MU(I.3.?) = (AMACH*B1 -.OMCGA) *1 
MU (2*3.2) = (AMACH*B?*OMFGA) *I*C 
MU (3.3*?) s OMEGA*0MFGA 
MU (1 .4*2) * (2.*0MEGA-.AMACH*B1 ) *1 
Ml) (2*4*2) (2. # 0MEGA*AMACH*B2)*I*C 
MUfl.S.2) = MU(1*4*2)*B1/I 
MU(?.5.2) = MU(?.4.2) *B2/I 
MUn.5.2) = B1»MU(3,2,2) 

Mu (1 .6*2) * C*92*MU(3.2.P> 

MU ( 2 * IS . 2 ) = CMPLX (2.0.0. 0)*MU(3. 1 . 2 ) 
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MU<?,6,2> = CMPLX (GAMMA /4«0»0,0) 

MU(3,6,2) * BETATMCMPLX (1.0*0.) *C) 

MU(2*H*2> = CFXP(I»B1»LENGTH) 

MU<3,8,?> = CFXP(I*R2*LEMGTH)*C 
MU<1*7,2) = RETAN* ( MU ( 2 , * , 2 ) ♦ MU ( 3, 5 » 2) ) 

MU<2*7*2> = OMEGA*0ETAI*(AMACH*<B1+C*B2> ♦OMEGA# (CMPL X ( 1*0«0.0) +C) ) 
MU <3, 7, 2) = OMEGA*BETAN* (AMACH* <B1* M U(2*B,2) *B2*MU < 3.8.2) ) ♦OMEGA# ( 
$MU(2»8,2) +Hlin ( 8*2) ) ) 

MU (1*9,2) = AETAC#BESSEL (MHAT+1,LHAT) 

MU<2»8,?> = I#MU(3g4.2> *(CMPLX(1 .-AMACH#AMaCH,0.0)*91-OMEGA*AMach) 
Munf8.2)=I#MIJ(3*4»?) #(CMPLX U.-AMACW*AMACH,0.0) #B2-0 MEGA#AMACw) #C 
MU(1*9»2) = -OMEGA* (AMACK#B1*0MEGA) 

MU(2* < 9,2) * -OMEGA* ( AMACH*B2 + 0MEGA) *C 

00 1 NX » 1 ♦ NDEX 

NBAR = NX - 1 

DO 3 LBAR = 1 , LOEX 

DO 1 L = 1 , LOEX 

Mu<L*21,2>=CMPLX(TNTFGRL(n,0.0»1.0,LHAT,L,L8AR,2,l> *0.0) 
MU(L,22,2)=CMPLX(XNTEGRL (I2*0.0»1.0.LHAT,L»LBAR,2,1> ,0.0) 
MU<L*23.2)aCMPLXUNTEGRL(I3*0o0»1.0.LHAT,L,LRAR,2*l) *0.0) 

MU (L ♦ 24 , 2) =CMPLX ( INTEGRL (I4,0.0»1.0»LHAT »L »LB AR , 2. 1 ) .0.0) 

MU(L*25*2) =CMPLX ( INTEGRL (I5,0®0.1.0,LMAT *L»LBAR,2» 1) *0.0) 
MU(L.26.2)=CMPLX ( INTEGRL ( 1 6 . 0 . 0 ♦ 1 . 0 , LHAT ,L , LBAR , 2 , 3) ,0.0) 
MU(L*27,2)=CMPLX( INTEGRL (I7*0.0«1.0,LHAT,L,LBAR,2.3) ,0.0) 
MU(L»28.2)=CMPLX(INTEGRL(rfl,0.0,1.0,LHAT.L,LBAR.2,3) ,0.0) 

MU (L » 29. 2) =CMPLX (INTEGRL (19 ,0.0, 1.0. LHAT, L. LBAR. 2, 3) .0.0) 

1 MU (L« 30 ♦ 2) =CMPLX (INTEGRL (IiO,0.0»1.0,LHAT,L»LBAR,2,3) ,0.0) 
MU(2*14,2> = CMPLX ( 0.0 , 0.0 ) 

MIJ(3.14,2) = CMPLX ( 0.0 , 0.0 ) 

DO 2 NY * 1 , NOEX 
N a NY - 1 

MU (3,9,2) = MU(3g2,2)*CMoLX(FL0AT(N) .0.0) 

MU ( 1 , 1 0 ♦ 2 ) = MU (3,9,2) *MU (3*9,2) 

MU (2* 10,2) a MU (3,4,2) *C-MPLX (FLOAT (N ) ,0.0) 

MU (3,10*2) a CMPLX ( (-1 . ) ##N ,0.0) 

MU ( 2 * 16,2) a MU (2, 15.2) *MU(3,9,2) 
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rtUO,15*2) 
MU ( 1 ,16*2) 
MU (1*11*2) 

$ 

MU ( 1 « 12,2) 
S 

MU ( ?» 12, 2) 

$ 

MU (3, 12,2) 
s 

MU ( 1 » 13*2 ) 
S 

MU (? * 1 3, 2 > 

* 

MU (3* 13*2) 
f 

MU ( 1 » 14,2) 
S 


Rl*MU(3,9,2)*I 

C*B?*MU<3*0.?)*I 

ZCR0SS1 (I*B1 ,N,N8AR,MU(3,4,?) *LENGTH) - 
ZCROSS 1 (I*B) *N*N9AR,MU(3»4»?) ,0,0) 
ZCROSSHI*B?,N,N0AR, MU 0,4,2) , LENGTH) - 
ZCR05S1 (I*B?,N,N3AR,mU( 3,4.2) ,0,0) 
ZCR0SS2(t*91 , N,NBAR *MU (3,4,2) ,LENGTH) - 
ZCR0SS2 ( I*R) ,N,NBAR,MU(3,4,2) ,0.0) 
ZCROSS? ( I *8? , N» NBAR *MU (3,4,?) »LENGTH) - 
ZCROSS? (I*H?,NeNBAR,MU(3»4»2> ,0.0) 
ZCR0SS1(I*P1,N,MBAR,MU(3,4.2),X2)- 
ZCROSS1 (I*H) ,N,NBAR,MU (3,4,2) , XI) 
ZCROSS 1 (I*R?,N,N3AR,MU(3,4,2),X2)- 
7CROSS1 (I*B^,N,NSAR,MU(3,4,2) ,X1) 
ZCROSS? ( I «R] ,N,N9AR,MU(3,4,?),X2)- 
ZCROSS? (1*81 ,N,NBAR,MU (3,4,2) , XI ) 
ZCROSS? (I*R?»N,N9AR ,MU (3,4*2) « X2 ) - 
ZCROSS2 (I*8?»N»N9AR, MU (3,4*2) ,X1> 


DO ? L = 1 , LOEX 

MU (l, 11, 2) = CMPLX(RESPRTM(MHAT*1,L) ,0.0) 

MU ( ? * 1 1 , 2 ) = MU(1,11,?)*CMPLX<RLMDA,0.0) 

Ml J (?* 14.2) = MU(?,14,?» ♦MU(L,NY,1>*(MU(2,1.2)*MU(L,22,2)* 

$ (MU (1,2,2)* (MU(3,l,?)*MU(3,lIt2)-MU(3.R,2)*Mu<2,12,?> ) 

$,MU (2,2,2)* (Mu (3,1,?) * M U (1,12,2) -MU (3,9,?) *M(J (3,12,2) ) 

$*MU ( 1,3,2) *( (MUO, 3,?) *MU( 1,10,?) >*MU( 3,n ,?) ♦MU(?,1*«2> 

S*MU(2»12»R) > ♦MU(2,3,?)*( (MU(3,3,2) * M U ( 1 ,10,?) )*MU( 1,12*2) ♦ 
SMU(2,16,?)*MU(3,12,?> > >♦< (RLM0A*MU(1 » 11 * ?) *MU (L , 25 , 2) -MH AT* 

I (MU (1,11*2) *MU (L » 23, ?) +PLMDA* MU(L,?4,2) > ♦CMPLX (?. *MHAT*MHAT , 0 . 0 ) 


$*MU(L,21 ,?) ) 

5* (MU ( 1 ,4,?)*MUO, 11 ,?) +MiH?,4,2) *MU ( 1 , 1 2,2) -MU ( 3 , 9 , ?) * (MU (2, 12,2) 
$*C*MU (3 , 1 ?» 2 ) ) > ♦ MU (2,1 0,2) *MU (L»2 ?,2) * 

$(MU(1,5,2)*MU(2,12,2)*MU(?,5,2)*MU(3,12,2)-MU(3,15,2)*MU(3,11,?)- 
SMU ( 1, 16,2 >*MU( It 12*2) ) )/?. 

$*MU(2,6,2) • < ( RLMOA-MU (1,11,2) *MU (L, 25 » 2 ) -MHAT* <MU< 1 , 11 «2) * 

S MU(L,23«?) ♦RLM[)A*MU (L,?4,2) > +CMPLX ( 2. *MHAT*MHAT , 0. 0 ) *MU <L*2 1 ,?) ) 
5* (MU (3,6,?) ♦MU (3,10,?) *Mtj( 1*7*2) ) +MU (L , 22, 2 ) * 



041 


$(MU(2.7.2>*MU(3*10*2>*MU(3*7*2) >- MU ( 1 * fl . 2) *PESSEL ( MHAT ♦ l *L ) * 

SBESSEL <l*L54P>*(Mu (2*8*2) *MU(3* 13*2) ♦ MU(3*B*2> *M(j ( 1 * 14*2) 
f*MU(l*9«2>*MU(l«13*2> *MU (2*9*2) *MU(2.l3*2> ) ) > 

2 MUn*l4«2) = MU(3*l4*2> *MU (L «NY * 1 > * (MU ( 2 ♦ 1 • ?> * M U(L .27 .2) * 

* (MU (1*2*2)* (MU(3*1*2)»MU(3» 11 *2) -MU (3*9.2) *M|.l ( 2, 1 2. 2) ) 

^♦MU(2*2*2)*(Mfj(3« I,2>*MU(1*12*2)-MU<3.9,2) *M(J(3» 12*2) ) 

S + MU ( 1 ♦ 3 * 2) * < ( MU (3*3*2) ♦MU (1*10*2) > *MU (3*11*2) ♦MU (2*16*2) 

$*MU(2* 12*2) > ♦MfJ (2*3*2) *< (MU (3*3*2) ♦MUM * 10*2) ) ««UU * 12*2) ♦ 
SMU(?,16*2)*MU<3, 12*2) ) ) ♦ ( (RLM0A*MU ( 1 * 11 * 2) *MU (L . 30 * 2) -MHA T* 
*(MU(1»11*2>*MU(L,28,2) ♦PI.MDA* MU(L*?9,2) )> 

$*(Ml)(l*4*?> *MU (3* 11*2) ♦MU (2*4*2) *MU ( 1 * 1 2 *2) -MU ( 3 * 9 . 2 > * < M <* ( 2 ♦ l 2 • 2 > 
5*C*MU(3* 12*2) > ) ♦ MU(2*in»2)*MU(L*27*2)* 

$ (MU (1*5 * 2) *MU (2*12*2) ♦MU (2*5*2) *MU (3*12*2) -MU (3*15*2) *MU ( 3 • II*?) — 
$MU(1 *16*2)*MU(1*12*2) ) >/?. 

5*MU (2*6,2) * { (RLMDA*MU( 1 *1 1 *2) *MU(L»30*2) -MH AT* ( MU (1*11*2)* 

$ MU (L«28»?) ♦RLMD4«MU(L*?Q*2) ) ) 

S*(Mll(3*6*?> ♦MlJ(3*10*?)*Mll(l *7*2) ) ♦MU(L,27*2>* 

* (MU (2*7,2) *MU (3* 1 0 *2) *MU( 3* 7, 2) )- MU < l * 8* 2> *BESSEUMHAT + 1 *L ) * 

SBESSEL (3*LBAR)*(MU(2*ft*2) »MU<3*13*2) *MU ( 3* B* 2) * M U ( 1 * 14*2) 
$+MU(i*9*2)*MU< 1*13*2)*MU(2*9,2)*MU(?*13*2) ) ) ) 

MU ( 3 *2 1 » ?) = A 1 ( 1 ) /FNORM (LBAR » 0 *NBAR , BESPPIM ( 1 *LR AR ) * LENGTH) 

MU (3*22*2) =A 1 (2)/FN0PM(LRAP*?,NeAR,8ESPRIM(3.LBAP) *LENGTH) 
MUC°0SS(LBAP*NX*1> = MU U , 1 * 2) *MU ( 3 *21 * 2) *MU ( 2 * 1 4 *2) 

3 MUC p 0SS(LF)AP,NX,2> = MU < 1 ♦ 1 * 2) *MU < 3 * 22 * 2) *MU ( 3 . 1 4 « 2) 

RETURN 

400 FORMAT** *50 ( 1 Ht ) /* WITH THE CROSS PRODUCT TERm$*/* *50{1H*>> 

END 


I 
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COMPLEX FUNCTION ZCROSSI (A.M,N,PI,X> 

OIMFNSION PI (?) 

COMPLEX A « ZB, ZD, ZF 

C ******** # *******************»#»»*<HHH»**«#*»*****#lHH»***<»** < HHMMHHHH M tJ Mt 
c ««« NOTE i IN THIS FUNCTION SUBPROGRAM PI IS DEFINED AS PI ( 3.14 ETC ) 
C *** OIVIOED BY LENGTH 

ZA (P.X) = SIN (P*X) /P/2, 

ZB (A »P, X ) s (A*STN(P*X) -P*COS(P*X) ) /<A*A+P*P) 

ZC (P»X) = (X*SIN (P*x)/P)/2. 

ZD ( A » X ) = { A* X~1 . ) /A/A 
ZE (P.X) = -COS (P*X) /P/?. 

ZF ( A «P , X ) = (A*COS(P*X> +P*StN(P*X) ) /(A*A+P*P) 

ZG (P.X) = SIN{P*X> **2/P/2, 

IF ( M ,EQ* N ) GO TO 1 

p A = (M-N) * PI S PR s (M+Nl * PI 

ZCROSSI = CEXP(A*X)*(ZA<PA,X) *ZA < PB » X > -A/2. « (ZB { A , PA , X ) /P A +ZB ( A , PR 
$,x)/P9) > $ return 

1 IF ( CABS ( A ) .LT. 1. E-lOf) ) GO TO 3 
IF ( N .EO. 6 > GO TO 2 
PA = 2,*N*PI 


ZCROSSI = CEXP(A*X)/?.*< CMPLX (1.3«0.0)/A*ZF(A.Pa«X) ) S RETURN 

2 ZCROSSI = CEXP(A*X)/A S RETURN 

3 IF ( N .NE. 0 ) GO TO 4 

ZCROSSI = CMPLX « X , 0.0 ) $ RETURN 

4 PA = 2.*N*PI 

ZCROSSI = CMPLX (ZC(PA,X> ,0.0) * RETURN 

ENTRY ZCROSS? 

IF ( M .EO. 6 > GO TO 7 
IF ( M «EO. N 1 GO TO 5 


PA = (M-N) * PI t PB r (M + N) * PI 

ZCROSSI = CEXPtA*X>*<ZE(PA,X>*2E<PB,X)*A/2.*<ZF<A.PA,X>/PA*ZA<PR, 
$)/PB> ) $ RETURN 

5 IF ( CA8S<A> .LT. i. E-IOn » GO TO G 


PA = 2 , *N*P I 

ZCROSSI = CEXP{A»X)*78(A, PA, XJ/CMPLX 12.0.0*0) 5 RETURN 


X 


6 ZCROSSI = ZG(N*PI,X) $ RETURN 

7 ZCROSSI = CMPLX ( 0.0 , 0.0 > $ RETURN 

END 
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complex function zntegrl (Nhat*n,a*piorl»x, type) 

INTEGER TYPE 
COMPLEX A 

P = NHAT * PIORL S 0 = N * PIORL 

GO TO ( 1 ,S) » TYPE 

1 IF(NHAT.EO.N.AND.N.EO.O) GO TO 2 

IF (NHAT .EQ.N. AND.N.EO. 0 • OR. NHAT. EQ. 0 ) GO TO 3 
IF ( NHA T • EO.N. AND .N.NP • 0 ) GO TO A 

ZNTEGRL - CEXP( A*X)*(CMPLX (SIN< (P-0> *X>- (P-O) ♦SINt (R*Q)*A)/(P*0) *0 
$.0) -A* ( (A*$tN ( (P-Q) *X) - (D-0)*COS (P-Q) *X) ) / ( A* A* { P-O* (P-Q) > ♦ ( A*S I N ( 
S(P>Q)*X)-(P*O)*C0S(( D *O>»X))/(A*A+(D*Q)#(P*.Q)>)) /CMPLX (2. 0 *0.0 ) 
RETURN 

2 ZNTEGRL = CEXP(A*X)/A S RETURN 

3 ZNTEGRL = CEXP ( A*X) / < A*A + P*P) * (A»C0S (P*X) +P*$IN (P*X) ) $ RETURN 

4 ZNTEGRL = CEXP ( A*X) / (A*A+4.*P*P) * ( ( A*C0S (P*X) *2 . *P*S I N ( P*X ) ) *COS (P 

S*X) *2.*P*P/A) $ RETURN 

5 IF ( N «EO* 6 ) GO TO 6 

IF (NHAT. EO.N. AND. NHAT.EO.OI GO TO 7 
IF(NHAT.EO.N/ANO.N.NE.O) GO TO 8 

ZNTEGRL = “CEXP (A*X) /CMPLX ( 2 . 0 ♦ 0 . 0) * (COS ( (G-P) *X ) / (Q-P) ♦C05( (Q*P)* 
$X)/(Q*P)-A«( (A*COS ( (O-P) *X) ♦ <Q-P)*SIN( (Q-P) *X)/( A*A + (Q-P) * (G-P) ) + A 
$*COS ( (G*P) *X) ♦ (Q + P) *S IN ( (Q + P) *X) ) /<A*A+ (Q+P) * (Q+P) ) ) > S RETURN 

6 ZNTEGRL = CMPLX (0. 0*0.0) <5 RETURN 

7 ZNTEGRL = CEXP ( A*X) * ( A*STN ( G*X) -Q*COS <Q#X ) > / ( A*A +Q*Q) S RETURN 

8 ZNTEGRL * CEXP(A*X)/CMPLX (2.0*0. 0)*(A»SIN(2.*P*X)-CMPLX(2.*P. 0.0)* 

SCOS(2.*P*X) )/(A*A*CMPLX(4.*P*P*0.0) > $ RETURN 

ENO 



. 044 


C 

c 

c 

c 

c 

c 

c 
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SU8R0UTINE EIGEN 
integer REDEFIN 

REAL LOWLIM* LENGTH* THETA(3)* INTEGRL 

REAL II* 12* 13* 14* 15* 16* 17* 18* 19* 110* Ill* 112* 113* 114* 
SI 1 5 

COMPLEX A* El* E?* E3» E4* 81C0N* 8?C0N« 63* ZINTGRL* Z<?4>. CCON. 
SFNOPM. PSI* WM91* WMR2* OMOD* BIMOO, R2M00* f/81M* WR?M. WMB ICON ♦ 
*WM8?CCN. pSISTR* PSICON, BETAI2, 0MEGA2 • OMFGAC* RET4TC* ANC* P(1B 
$)*ZINT(41, REVAL(5), ZNTFGRL. SLOPE* RETASTR 
COMPLEX 81 * 82 * BETAN* 8ETAC ♦ 8ETAI* BET AWT , OMEGA* C* AK 1 * 

1 AK? ♦ PSI?* MU* BETflTS. Cl?* I* MUSNO 
COMMON / BLKfl /HI * R? * AMACH, GAMMA* OMEGA. ALENGTH, c* AK 1 
1* AK2, BETAN, RLMOA* BETMt BETAC* BFTAWX* LhaT* MHAT* NHAT, 
1PSI?*X1*X? 

COMMON / RLKR / MU(3*30»?) * 8ETAIS(2) *1 * MATRIX 
COMMON / RLKC / LENGTH * MUSNO ( 30 * 3*2) * Al(2) 

COMMON / BLOCK / C 12 ( 30 ♦ ^0) * I SET .LDEX ,NDEX .ERROR I .ERR0R2 
COMMON / ROOT / 8ETAT2 ♦ RETASTR . SLOPE 

EQUIVALENCE <R ( 1 ) *M(J ( 1 * 1 * 2> > * ( Z <1 ) *MU ( 1 * 7 * 2) > , (Z INT ( 1 > *MU ( 1 , 15 
$.2M * (PEVAL(l) *MU<2*16,2)> 

EXTERNAL II* 12* 13. 14, IS* 16* 17, 18* 19, Ilf), III* 112, 113* 
$114* IIS 

EXTFRNAL Rl* R2» R3, R4* R5* R6* R7* R 8* R9 
DATA PI / 3.1M59265358970 / 

ZINTGPL < A *P* UPL IM.L0WL1M ) = <CEXP < A*UPL I M) * < A*COS ( P*UPL I M ) + P*S IN ( 
fP*UPLIM) )-CEXP(A*LOWLIM)*(A*COS(P*LOWLIM)*P*SIN(P*LOWLIM) ) ) / ( A*A*P 
S*P) 


000000 » 0 *ft# 0»0400 

00000000 #** 0 * 0**0 

00000400»«»0ft0#00 

00440000*#***0#*0 

04040000*»»*»0**4 


0400000404044 

NOTE. IN THIS ROUTINE THE INTERACTION ****»»**«#**» 
INDEX <N> HAS BEEN REPLACED RY THE ********»**#» 

ALPHA CHARACTER B* SO AS NOT TO 8F *-*«•**#**»**** 

CONFUSED WITH THE MATRIX INDEX IN) * ***»»**«#**«* 

00000»0»4»00» 


00«04040»«*0»000040400040440040»«»«0440404»«4444»4#04000440»0«00000400 


81C0N = CON JG ( B 1 ) 
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8 1 MOD = CABS(81)**2 
8?C0N » CONJG (82) 

B?MOD = CAaS(B2>*»2 
B3 = 81 ■> 8? 

CCON = CON JG < C > 

CMOn = CABS i C ) **? 

const = nhat * pi / length 

El = I* <2o®fll“BlC0N) 

E? = I*(2o*B?~82CONji 
E3 = I*fll 
E4 s 1*82 

0MEGA2 = CM°LK (0o9*0»0) 

PIORL = PI / LENGTH 

PSISTR = PSI? * ENORM (LhaT*MHAT, NHAT, RLMDA, LENGTH) **2 
PSICON = C0NJG(C5QRT(PSISTR) ) 

W01M = OMEGA ° AMACH * 01 CON 

WB2M = OMEGA » AMACH * BPCON 

WM8 1 * OMEGA « AMACH * 8) 

WMB2 = OMEGA * AMACH * HP 

WMBlCON = OMEGA « AMACH * BlCON 
WMBpCON = OMEGA ® AMACH « B2C0N 
WRITE (fe«602> 

C ***************** 

C ft**************** THESE ARE THE PHIU) ** 3 PRODUCTS »»«««»«»»»«»««« 
C »***»* »»**#«»*»*«• **»**»<>»«*»»**« 

C «*****»•»««**«*«•****■»■*■»•»*•** »**•»*#•»*•& ««»««» a «««««« »***#**•»*. 

C *** 

C *** THESE ARE THE THETA-INTEGRALS 
C *** 

THETA ( 1 ) = 75 * PI 

THETA (2) a PI / 4 0 
THETA (3) = PI 

C *#» 

c *** these ARE the R-INTEGRALS 

C *** 
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Z(l) = CMPLX ( INTEGRL {Rl.O*0.1*0.LHAT»l*l*2.2) *0.0) 

Z f 2) = CMPLX<INTEGRL<R2.0.0.1.0.LHAT.1.1.2.2> .0.0) 

Z ( 3) = CMPLx<INTEGRL<R3.0.0.1.0.LHAT.I»1.2.2) .0.0) 

Z<4) * CMPLX(INTEGRL(R4.0.0.l.O»LHAT.l *1*2.2) .0.0) 

Z (5) = CMPLX (INTEGRL (R5.0, 0.1*0. LHAT»1.1*2*?) » 0 . 0 ) 

Z (6) * CMPLX(tNTEGRL(R6,o,0.l.0.LHAT.l. I.2.?) .0.0) 

Z ( 7 ) = CMPLX(!NrEGRL<R7.O.0»l.0.LHAT,l.l.2.2) . 0 . 0 ) 

Z (8 ) = CMPLX(INTEGRL(Ra,0.0.1.0.LHAT.l.l,2.2) .0.0) 

R < 1 ) = RLMOA*RLMOA*Z (1 >-P.«’MHAT*RLMr)A*Z (?) ♦MHAT*MHAT«7 (3) 

R (2) = Z(3> 

R(3> = Z (4) 

R (4> = RLMOA*RLMGA* (-RLMOA*RLMDA*Z ( 1 > ♦ 2 ,»MH A T»RLMQA»Z (2) -MHAT*MWflT 
S*Z <3>-RLMDA*7(5) ♦MHAT*(3* *MHAT)«Z<6> ) -MHAT*MHAT* ( 3. ♦P.OMHAT ) «Z ( 7 ) * 
SRLMOA 

R (5 ) = RLMD4*RLMDA*Z<6>-?.*RLMDA*MHAT»Z(7)*MHAT*MHAT*Z(R) 

R (6 > - RLM0A*Z(7)-MHAT*Z(fl) 

R (7) s Z(fl) 

. R ( 8 > = CMPLX (R9 { l * » RLMOA .FLOAT (LHAT ) .FLOAT (MHAT ) ) * 0 • 0) 

REVAL (1) a CMPLX (RLMDA*RLM0A*R1 (1.0.RLMDA,0.0»0.0)-2.*MHAT*RLMDA» 
SR2(1 .0. RLMOA .0.0 .0.0) ♦MHAT*MHAT«R3 (1 .0 . RLMOA ♦ 0 . 0 . 0. 0 ) .0.0) 

REVAL(2) = CMPLX(R3(1.0. RLMOA. 0.0. 0.0) .0.0) 

REVAL(3) a CMPLX {R4 { 1 .0 *°LMDA .0 • 0.0.0) *0.0) 

REV4U4) = CMPLX(R9(I .0. RLMOA. 0.0.0. 0) .0.0) 

C *«* ' 

c these are the z-integrals 

c 

ZINT(l) * ZINTGRL ( El .CONST. LENGTH. 0.0> 

Z1NT(2> * ZINTGRL C E2.C0NST. LENGTH. 0.0) 

ZINT ( 3 ) = ZINTGRL ( E3.CONST.LENGTH.0.0) 

ZINT<4> s ZINTGRL ( E4. CONST. LENGTH.0.0) 

Z(l) = -<WMRl*WM01*(7INT(l)*ZINT(3)*(C+CCON)*CMOD*ZINT(4)). 

% WM82*WM82* (ZINT <3> * ZINT ( 4 ) * (C+CCON) ♦CMOO*ZINT (2) ) ) 

Z (2) = -(WM81*WMR1* (HIM00*ZINT (1 > ♦ZINT (3? * <C*82*RlCQN.8i*CC0N*R2C0 
SN» *CM00*8?M0D*ZINT<4) ) + 

l WMH?*VMB2* (B2MOn*ZINT (3> ♦ZINT(4) * (C»b2*RlC0N + Bl*CC0N*R2C0 

SN) ♦CMOO^BPMOO^Z I NT (?) ) ) 
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Zd) = -(WMR1*WM81*(ZINTU>*2.*C*ZINT(3)+C«C*ZINT(4>) + 

<5 CC0N*WMR2*WMB2*(ZINT (3) *2«*C*2INT (4 > *C*C W Z INT (2) ) > 

Z (4 ) = WMRl*WM81*<Bl*RI»7lNT f 1> 42.*B1*B2*C*ZINT<3> ♦C*02«C*82 *Z INT i 

S4) ) 4 

S CC0N*WM8?*WMB2*(Pl!*Bl*ZiNT(3) ♦2.*R1*82*C*ZINT (4) 4C*B2*C*02* 

SZINT (2) ) 

Z (5) = ZlNTt JM(2»*CoCCOfr!)*ZIMT(3)+C«(2.*CCON4C)*7lNT(4)4C*CMOr>*Z] 
?NT(?) 

Z (6 ) = 81C0NMB1*ZINTU>oC*B 3*ZINT<3>*C*C*R?*ZINT<4>>4 
$ CC0N*fi?C0N*(Bl»ZINT(3> ♦C*83*Z1NT (4) *C*C*H2*ZINT (2) > 

Z (7) - -<Sl*Bl*ZINTO) *B1*(2.*C*H24RI*CC0N) *7 INT (3) + 82* (C*C*82«-2. « 
5B1*CM00) *7lNT (4) ♦C»CM00*R2«82*ZINT (21 ) 

Z (8) - -(RI»B1*81M00*ZINT(1> ♦<C*82*B3*BIMOO4CCON*82CON*0I*B1*B1 ) * 
SZINT(3> ♦ (C’»C*B 2 *B 2 * 9 ?*BlCON 4 CMOO'* 8 l*B 3 * 82 MOr))*ZINT( 4 ) ♦C*CM 0 D*B?*R 2 
$*32MOO*ZINT(2>) 

Z (9) = 81 HOD*B 1 M0D*Z INT (1>+(2»*C*B1 MQ0*B1 C0N*82+CC0N*Bl*8l *B2C0N*8 
S2C0N)*ZINT (3) ♦ (C*C*B2*B2*filCOM*BlCON*2.»CHOD*B2MOD*B2CON*Bl ) »Z*NT ( 
$4> ♦CM9O*CMOO*02MOO*B2MOD»ZINT (2) 

Z(1G> = -( (0MEGA>AMACH*Bi/2. )*ZINT(3)+C»(0MEGA44MACH*B2/2.)*ZINT(4 
t) ) 

c *** 

C »«* REDEFINE THE INTEGRATION LIMITS FROM XI TO X2 
C *** 

IF ( CABS (RET AC) .LE* 1 <> E-l 0 »0R. X2“X1 .LF.l.E-10) GO TO 2 

ASSIGN 2 TO REDE^IN 

ZINTU) a ZINTGRL (El * CONST » X2*X1 > 

ZINT ( 2 ) ~ ZINTGRL (E2* CONST * X2 * X 1 ) 

ZINT (3) = ZINTGRL (E 3* CONST « X2* XI ) 

ZINT (4 > = ZINTGRL (E4« CONST <» X2* XI ) 

1 Z(ll> = I«(WMRI*WM81*WMB1CON#ZINT(1>4WMR1*(2.*C#WMBICON*NMB24CCON* 
SWMB1*1«M82C0N>*ZINT(3) ♦WMR2 *C w (2.*WMR1*WMB2C0N*C*wM82*WMB ICON) *?INT 
$ (4) +C*C*WM8?*WM8?»WMR2C0N*ZINT (2) ) 

Z ( 1 ?> = -I*(WMBICON*ZINTU) * (2.*C*WMR1C0N+CC0N*WMB2C0N)*ZINT(3) +C* 
S(C*WMB1CON*2.*CCON*WMR2CON>»ZINT(4) ♦C*CMOO*WMB2CON*ZINT(2) ) 

Z < 1 3 > = I*(Rl*8I*WMRlC0N*ZlNT(n4 ( 2 .«B 2 »C*WMR 1 CON* 81 *CCON*WMR 2 CON) 
$*ZINT (3)*81+B2*(WMeiC0N*C®C*82*WMB2C0N»CC0N*2.*Rl*C)*ZINT (4) *C*WMB 
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S2C0N*CM0D*e?»B?*ZlNT(?> ) 

Z(l4> = ZlMT (3) ♦C*ZINT (41 

Z (22> * B 1*7 TNT ( 1 > ♦ ( C»H3*CC0N»8 l ) *ZTNT (3) ♦ (CM00*«3*C*C»B?> *ZINT (4) 
$*C*CM0D*B?*ZINT (?) 

Z(21> = -T*(R]CON»ZINT(l ) ♦ (CMPLX (2. 0.0. 0) *C*R1C0N ♦CC0N*B 

S2C0N )*ZINT ( 3 ) ♦ (CMPLX (2.0*0. 0> *CM0D»82C0N.C*C *B1C0N) 

t*Z I NT ( 4 ) *C *CM00*fl?C0N«7INIT (2) ) 

Z(24) = B!*B1M0D*ZINT ( 1 > ♦ (CMPLX (?,0. 0.0) *C*R?«8lM00*CC0N*82C0N*81 » 
$81)*ZINT<3) ♦ <C*C*8?*8?*B1C0N+CMPLX (?.0»0.0)*CMOO*B2MOO*R1 >*ZINT(4) 
?*C*«2*CMOn # B2MOO*ZlNT (2) 

GO TO REOEFIN 

C **» 

C *** EVALUATE THE VOLUME INTEGRALS 
C *** 

C •»** TERMS 4 AND 5 
C *** 

2 BETAI2 = ((GAMMA-1. )*(2.*(Z(1)*(THETA(1)*R(1)-»THFTA(?)*MHAT*MHAT* 
SR(2>> ♦THETA(1>*R(3)*Z(2) ) * 2 ( 3) * (THETA (1 ) *R ( l ) ♦MHAT*MHAT*THET A (? ) * 
SR (2) > * THETA (1)*R(3>*Z(4) > 

C *** 

C *«* TERM 6 

c *■»* 

*♦ (2»*Z (5)* (THETA ( 1 >*R (4) ♦MHAT*MH AT*THET A ( 2) * ( ? . *R ( 5) -R (6> -MHAT* 

$MHAT*R(7) ) )*Z(f.)*(THETA(l)*(3. ) *P ( 1 ) *3 . *MH A T*MHAT*THE T A ( 

S2)*R(2) ) *Z (7) * ( THET A(1)*o<1) *MHAT*MHAT*THET A (2) *9 (2) ) -THETA < 1 > *R (3 
S)«(Z(8)*Z<9) ) ) 

C *** 

C TERM 7 

C *** 

*♦ (Z (S>*(THETA(D*P(4) +MHAT*MHAT*THETA(2)*(2.»R(5)-R(6)-MHAT*MHA 
$T*R<7> ) ) ♦Z(A)»MHAT*MHAT*THETA(2)*R(?) *Z (6) * (THETA (1 ) *R ( l ) ♦MHAT*MHA 
ST*THETA(2)*R(?) >*THETA(1)*R(1)*Z(7) + Z (0) *THETA ( I > *R<3) ) > /CMPLX <16, 
S.f).n)/PSISTR/PSICON 

SLOPE = THETA (3)*R(8)*Z(1Q)/CS0RT(PSISTR)/I 
BETA 12 = - BETAI? 

C *»* 
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C *** EVALUATE THE SURFACE INTEGRALS 
C *** 

C *** FIRST DO THE LINER INTEGRALS 
C *** 

IF < CABS(8ETAC> ,LE.l.E-10.OR.X2-Xl.LE.l.E-10> GO TO 5 
BETAl? .a BETA 12 ♦ GAMMA*BETAC/CMPLX (2. 0*0.0) *( (-THETA ( 1) »REVAL*Z < 1 
SIM ( THETA (1) *REV AL ( 1) ♦MHAT*MHAT*THETA <2) *REVAL (2) ) *Z ( 1 2) ♦ THETA (I >* 
SREVAL (3>*Z(13)+ 

S(GAHMA*l.)/CMPLX (rt. 0.0.0) * (THETA <1)*REVAL (3> * ( I*OMEGA*OMEGA*OM£G4* 
$Z(5> ♦OMEGA *OMFGA* AMACH* (CMPLX < 2. 0 « 0 . 0 > *Z < 22 ) -Z (23> > ♦ I "OMEGA* AM ACH* 
$AMACH*(CMPLX(2.0*0.0)*Z(G)~Z(7) ) +AMACH* AMACH*AMACH*7 <?4) > > )/PSISTR 
S/PSICON) 

SLOPE = SLOPE ♦ GAMMA*BETAC/CMPLX(2.0*0.0)*(THErA(3)"REVAL(4)*?(14 
$)/CONjG (PS I CON) ) 

C *** 

C **»EVALUATE THE INJECTOR INTEGRALS 
C *** 

5 ASSIGN 3 TO REDEFIN 

ZINTU) = ZINT (2) - ZINT (3) = ZINT(4> = CMPLX ( 1.0 . 0.0 > 

GO TO 1 

3 GET A 1 2 = BETA!? ♦ GAMMA*RETAl/CMPLX (2.0 . 0 . 0 ) * < ( ‘THETA ( 1 ) *R ( 3) *Z < 11 
%) .(THETA (I)*R(1) ♦MHAT*MhaT*THETA(2)*R(2) )*Z(12) +THETA ( 1 ) *R<3) "Z ( 13 
S>* 

$ (GAMMA ♦!,) /CMPLX ( 8 . 0 . 0 . 0 > * ( THETA (1 ) *R <3> * ( I *OMEGA*OMEGA»OMEGA*z (5 ) 
$.OMpGA*OMEGA*AMACH* (CMPLX (2. 0*0.0 )*Z (2?) -Z (23) ) ♦ I "OMEGA* AM ACM* AM AC 
?H* (CMPLX (?.0.0.0>*Z( 6) -7 (7) ) ♦AMACH*AMACH*AMACH»7 (24) ) ) )/PSISTP/PS 
$ ICON) 

SLOPE = SLOPE ♦ GAMMA*HETAI*CMPLX(2.0.0.0)*(THETA(3)*R(8)*Z(14>/CO 
$NJG (PSICON) ) 

C *** 

C *** EVALUATE THE NOZZLE INTEGRALS 
C *«* 

ASSIGN 4 TO PEDEEIN 

ZINTU) = CEXP(Et*LENGTH)*CMPLX( (-1,)**NHAT.0.0) 

Z INT (2 ) = CEXP(E?*LF_NGTh>*CMPLX ( < - 1 . ) **NHAT . 0 . 0 > 

ZINT (3) = CEXP(E3*LENGTH)*CMPLX ( (-1. )**NHAT.O.O) 
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ZINT(A) = CEXP(E*»LENGTH1 *CMPLX< (-l.)#*NHAT,0.0) 

GO TO 1 

U BET A 1 2 = BETAI 2 ♦ GAMMA*RET AN/CMPLX < 2. C * 0 . 0 > * (( -THETA < 1 ) *R ( 3) *7 ( \\ 
«)♦ (THETA U)*R(1) +MHAT*MHftT»THETA (?) »P (?) ) »7 C 12) ♦ THETA (1 ) *P (3) *7 ( IT 
S ) ♦ 

* ( GAMMA* 1 ,) /CMPLX <fl. 0*0.0) * (THETA (1 ) *R (3) * ( I*0MEGA*0MEGA*0MFG4 ° 7 (5> 
?*OMEGA*OMEC)A*AMACH* (CMPLK (2.0*0.9>*7(?2>-Z<23>) ♦ I *OMEGA*AMACH° A mac 
?H* (CMPLX (?. 0*0*0) *Z< 6) -7(7) ) ♦ 4MACH*AMACH*AM 4CH*7 (?A) ) ) ) /PS l ST^/PS 
SICON) 

SLOPE = SLOPE ♦ GAMMA*BETAN*CMPLX(2.0«0.0)*(THETA (3) *R (0) *Z < 14) /CQ 
SNJG(PSICOM) ) 

Q **»*»»«*#»*#»***«• 

C **»*»»*»*»***«#»* THESE ARE THE PHI (2) , PHI(1) PRODUCTS *****##»#*« 
£ «*#»*»«»«»« 

00 15 L = 1 * LOEX 

C *** 

C *** THE THETA INTEGRALS REMAIN THE SAME AS THOSE DEFINED ABOVE 
C *** 

C **» THESE ARE THE R-INTEGPALS 
0 *** 

H ( 9 ) = CMPL X ( RLMDA* TNTEGRL ( 1 1 1 * 0 - 0 ♦ 1 . 0 .LHAT *LH AT ,L *MHA T ♦ 1 . 1 > -«H A T 
**INTEGRL (T 1?. 0.0*1. 0«LHAT*LHAT*L*MHAT*I*1) .0.0) 

R(lQ) = CMPLX (RLMDA* INTFGRL (113*0.0*1.0 *LH AT *LHA T»L«MHAT*1 * 3 > -MHA 
$T* INTEGRL (114.0.0*1 .0 *LH AT * LH AT * L * MHAT *1.3) *0.0) 

R ( 1 1 ) = CMPLX ( INTEGRL (11.0.0*1. 0 * LHAT *LHAT*L*MHAT*I»l) *0.0) 

R ( 1?) = CMPLX(INTEGPL(I6*0.0.1.0*LHAT*LHAT.L*MHAT*1.3) *0*0) 

R ( 1 3> = CMPLX (INTEGRL(I?.0. 0*1.0* LHAT, LHAT*L*MHAT*1*1) *0.0) 

R(l4) = CMPLX(INTEGRL(I7.0.0*1.0*LHAT*LHAT*L.MHAT*1.3> .0.0) 

R ( 1 S) = CMPLX (HESPRIM( l,D*( INTEGRL (1 15. 0.0 *1.0 *LHAT.LhTT.L.MHaT* 1 
$ » 1 ) -BESPR IM ( I ,L ) * INTEGRL (I?«0. 0*1.0* LHAT *LHAT »L*MHAT*1 *1 ) ) *0.0) 
R(10> = CMPLX (-BESPR I M (3*L) * (INTEGRL (115.0.0*1 • Q »LHAT .LHAT . L *MHAT* 
$1*3) ♦BESPPIMn,L>* INTEGRL (12*0.0*1.0* LHAT *LHAT*L*MHAT *1*3) ) *6.*INT 
SEGRL ( 10*0.0* 1.0*LHAT.LHAT,L*MHAT*1*3) *0.0) 

R( 17) = CMPLX (-BESPRTM ( 1 .L) *1 NTEGPL (T15»0» 0*1.0* LHAT «LHAT «L,MHaT + l 
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? i 1 ) .0.0) 


R(1S> = CMPL X ( BESPRIM<3.L>*INTEGRL(I15.0.0,1.0,LHAT.LHAT,L.MHAT* 
S.3>-2.0*INTEGRL(I6»0.0, I .0,LHAT,LHAT,L,MH4T*1.3) .0*0) 

REVAL(l) = CMPLX(RLMDA» Til (l.O.RLMOA.RLMDA.BESPPIMfl .L) )-MHAT»Il 
$(l.r .RLM0A.RLM0A.8ESPRIMn.L> ) .0.0) 

REVAL(?) = CMPLX(RLMOA* T1 3 (1 . 0 . RLMOA , RLMOA .BESPRI M (1 .L ) > -Mn A T* 1 1 
5(1.0. RLMO A .RLMOA , 8ESPR IM { 3 , L ) ) .0.0) 

REVAL(l) = CMPLX(I6(1.0.PLMDA.RLMDA,RESPRIM<3.L) ) .0.0) 

REVAL(4> = CMPLX (12(1 .O.PLMDA.RLHOA.BESPRIM (1 .L) ) .0. 0) 

REVAL(S) = CMPLX(I7<1.0»PLM0A, RLMOA, 8ESPRIM(3.L) > ,0.0) 

DO 15 NX = l ♦ NOEX 


N = NX 
ZINT(I) 
5 

ZINT (?) 
5 

ZTNT (3) 

s 

ZINT (A ) 
5 


1 

ZNTEGRL (NHAT .N, - 1*0 ICON.PTORL. LENGTH, 1 ) - 
ZNTFGRL (NHAT.N,-I»81CON,PIORL,0.» 1 ) 
ZNTEGRL (NHAT .N .-I*B2C0N , PIORL .LENGTH, i ) - 
ZNTEGRL (NHAT »N.-I*B2CON , PIORL , 0. « 1 ) 
ZNTEGRL (NHAT .N«- 1*0 ICON .PIORL .LENGTH, P ) - 
ZNTEGRL (NH AT »N. -I *B1C0N. PIORL. 0., 2) 
ZNTEGRL (NHAT » N. - I *B2CON» PIORL, LENGTH.?) - 
ZNTEGRL (NH AT. N,-I*B2CON,PIOPL.O., 2) 


C *** 

C *** THESE APE THE Z VOLUME INTEGRALS 
C **« 


Z < 1 ^ ) = I*(NP1M*ZINT (1>+0C0N*WB2M*ZINT(2) )-AMACH*PIORL*N*(ZINT (3) , 
5CC0N*ZINT (4) ) 

Z(16> = (WB1M*B1C0N*ZINT(3) ♦CC0N*B2C0N*WB?M*ZINT (4) > -AMACH*PIORL 
S*N*(BIC0N*ZINT(1) ♦CC0N*B?C0N*ZINT(2) >*I 
Z(17) = -I*OMEGA* 2.0 * (WMB1C0N*WMR1C0N*ZINT ( 1) ♦CCON*RMB? 

$CON*WM82CON*ZINT (2) ) «-AMACH*PIORL*N* ( WMB1C0N*WM91C0N*ZINT (3) +CCON*h 
5MB2CON*WMB?C0N*ZINT (4) ) 

Z ( IB) = -I* ( (OMEGA- AM ACH*B1 ) *ZINT (1 ) ♦ (OMEGA-AMACH*82> *ZINT (2) *CCON 
5) 


C *** 

c *** these ARE THE Z SURFACE INTEGRALS 
C *** 


IF ( CABS(BETAC) .LE.1.E-10.0R.X2-X1.LE.1.E-10) GO TO 12 
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ASSIGN 

ZINTtl) 

S 

Z INT (2 > 
S 

ZINT (3) 

$ 

ZINT (4) 

$ 

11 ZUQ> = 
$ 


12 TO REOEFIN 

* ZNTEGRL(NHAT»N.-I*BlCON,PIORL.X2.1>- 
ZNTEGRL(NHAT.N.-I*B1C0N,PI0RL.X1,1> 
s ZNTEGRU (NHAT*N,-I*B2C0N.PI0RL. X2.1 > - 
ZNTEGRL(NHAT*N,-I*B2C0N*PT0RL*X1 *1 ) 

=• Z NT EG PL (NHAT*N,-I*B1C0N,P10RL*X2*2>- 
ZMTEGRL (NHAT.N,-I«B1C0N,PI0RL.X1 .2) 

= ZNTEGRL<NHAT»N.-I*R2CON.PIORL.X2,2>- 
ZNTEGRL (NHAT,N.-I*62C0Nf PI0RL.X1 *2) 
OMEGft*(wMBlCON*ZINT(l) *CC0N»WMB2C0N*Z IN T (2) ) ♦ 
I*N*AMACH«PIORL/p,« (W'181C0N*ZINT <3> ♦ CCON*wMR?CON*ZINT <4> ) 


Z ( 2 ^ ) = ZINTU>*CC0N*7INT(2> 

Z(2D = “I* (R ICON* ZINT ( 1 ) ♦CCON*B2CON*Zl NT (2M 
GO TO REOEFIN 


c *** 

C **« EVALUATE THE VOLUME INTEGRALS 
C *** 

12 BET A 1 2 = RETAI2-(MUSND(NX.L.I)*THETA(3>*( (R(o)*ZUS)-N*PtOPL*Rni> 
S*7(l6>)/ 2.0 ♦ (GAMMA-1. >/4.0 * <R ( 13) *Z < 1 7 ) ♦ 7 . ( 1 » ) 

$* (R { IS ) ♦R(17)-N*N*PI0PL*pI0RL*R<13) ) > ) ♦MUSNO (NX .L .2 ) « ( THETA ( 2) * ( (R 
S(10>*Z(15>~N*PIORL*R(14>*Z(16> »/ 2.0 ♦ (GAMMA-1 .) /4. 

$ *(R(|4)*Z(17)»Z(18)*(R(16) ♦R(lfl)-4.*R(12)-N*N*PI0RL»PI0RL»R(1 

$4 ) ) ) ) ♦ THET A (3) /?• *MHAT *R (1 2 > *Z ( IS >> > /PSICON 

C *** 

C *»* EVALUATE THE SURFACE INTEGRALS 
C *** 

C *** FIRST DO THE LINER INTEGRALS 
C *** 

IF ( CAHS(BETAC) .LE.1.E-10.0R.X2-X1.LE.1.E-10) GO TO 13 
BETA 1 2 - BETA 1 2 * BETAC* GAMMA/2.0 /PSICON* (MUSNO (NX »L . 1 ) 

$*THETA(3)*(REVAL(4)*Z(1Q)-(REVAL(1)*Z(20)-N*pIORL*REVAL(4>*Z(2!) ) 

$ 2.0 )+MUSN0(NX.L.2)*(THETA(2)*(REVAL(S)*Z(1R)-(REVAL(2) 

SZ(2n)-N*PlORL*REVAL(S)*Z(?l) )/ 2.0 )“ MHAT»THETA (3) / 

S2.0 *REVAL<3>*Z<20m 

C *** 

C *«* EVALUATE THE INJECTQP INTEGRALS 


* \ 



\ 
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C »*• 

13 ASSIGN 14 TO REDEFIN 

ZINT(l) s ZINT (2) - ZINT (3) = ZINT(A) a CMPLX ( 1 . 0 . 0 . 0) 

GO TO 11 

14 BET A 12 = BETA!? ♦ BETAI* GAMMA/2oO /PSICON* (MUSNO < NX «L • 1 ) 

S*THFTA(3>*(P(13>*Z(19>-(P(9>*Z(20)-N*PIORL*R(13>*Z(21> >/ ?.0 

$ ) +MUSNO <NX*L*2)* (THETA (2) *(R ( 14) *Z (1 9 ) “ (R ( 1 0 ) *Z ( 20 ) -N*PIORL* p ( 

$14>*Z<21> ) >- MhAT*ThETA (3J/2.0 *R(12)*Z(20) ) ) 

C *** 

C *** EVALUATE THE nozzle INTEGRALS 

c *** 

ASSIGN is TO REDEFIN 

ZINTd) = CEXP ( ”1*0 1 CON) * ( -1 . > *»( NH AT *N) 

ZINT (2) - CFXP(-I*B2C0N)* (-1 . ) ** ( NHAT+N) 

ZINT CD = CEXP(-I*81C0N) * ((“1* ) **NHAT* U (-1 . ) **N> > 

ZINT (4) = CEXP (-I*B2C0N) » ( < -1 . ) **NH AT* ( 1 •- ( -1 . ) **N) ) 

• GO TO 11 

15 6ETAI2 = RETAI2 ♦ BETAN* GAMMA/?. 0 /PSI CON* (MUSNO (NX «L * 1 ) 

S*THETA<1)*<R(13)*Z(191-<P<9>*Z(ZO>-N*PlOPt*R<13>*Z(21 ) »/ ?.0 

$ ) ♦ MUSNO (NX ^*2)* <THETA(2>*(R<14)*Z(19)-(R<10)*Z<?n>-N*PlORL*P< 

$14)*Z(21) ) )- MHAT*THETA(3)/2.0 *R ( 1 2) *7 (20 ) ) ) 

flETASTR = 8ETAI2*EPSIZN(LHAT»MHAT*NHAT»LENGTH.RLM0A)*PSI (0 >/FNORm ( 
SLHAT »MHAT .NHAT »RLMDA » LENGTH) /AK 1 
SLOPE = (SL0PE*BETASTR/BFTAI2) *1 
WRITE (6*700) SLOPE 

CALL SOLSION ( OMEGA * BFTAI ♦ GAMMA ♦ AMACH * °I ) 

*02 FORMAT (*©*///> 

700 FORMAT ( *0 THF SLOPE IS *PG21.14> 

RETURN 

END 


f 
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SUBPOUT INF NT AU ( AN ♦ OMEGA ) 

REAL N , NVEC 

COMMON / NORM / N(3) . TAU<3> » NVEC ♦ TAUVFC * PRINT 
OIMENSION AN ( 2 > 

PIE2 = 2. * 3. 14159245350P79323846 
ANR = AN ( l > 

AM = AN < ? > 

N s ( ANR**? ♦ AM**? > / 2. / ANR 
IF { ANI .LT. 0.0 ) F, 6 

5 TAU = (PIFP-ACOS < 1 .-ANR/N) ) /OMEGA 
GO TO 7 

6 TAU = ACOS ( l.-ANR / N > / OMEGA 

7 IF ( PRINT .EG. 10HNO PRINT ) RETURN 
WRI TE (6 » 1 00 ) Nil) . TAU ( 1 ) 

100 FORMAT(*0 THF INTERACTION INDEX IS *G2l.I4»?3X* THE SENSITIVE T I MF 
S LAG IS * G21.14) 

RETURN 

entry NORMAL 

OFT OF A s TAU(1)*TAU(1)*(TAU(2)-TAU(3>) ♦ TAU (?) *TAtJ (2) * ( TAU ( 3 > - T 
SAU ( 1 ) ) ♦ TAU(3)*TAU(3)*(TAUU)-TAU(2) > 

ANR = (N(1)*(TAU<2)-TAU(3) > ♦ N (2) * ( TAU < 3 > -TAU ( 1 ) > ♦ N<3>*<TAUU>- 
STAU (2 ) ) ) / OET OF A 

AM = (TAU< l)*TAU(n*(N(?)-N(3) ) ♦ TAU (2) *TAU (2) * <N ( 3) -N ( 1 ) > ♦ TAU 
S(3)*TAU(3)*(N(1) -N (?) ) > / DET OF A 
ANR = 2.*ANR*TAU (2> *ANI * PIE2 = SQPT ( ANR* ANR* 1 . ) 

TAUVEC = ANR / PIE2 * NVEC *» -1. / PIF2 
IF ( TAU .LT, TA'J (2) > GO TO 1 
NVEC - “NVEC $ TAUVEC = - TAUVEC 

l WRITE (6. 600 > NVEC * TAUVFC 
GOO FORMAT <*0 NVFC = *G?1.1A* TAUVEC = *G?1.14) 

RETURN 

END 
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SUBROUTINE SOLSION ( OMEGA * BETA I . GAMMA . AMACH , PI ) 

REAL NVEC * N? * LOwLIM 

COMPLEX OMEGA . 8ETAI • BET A 1 2* OMEGA? ♦ OMEGAC , BET A IC . ANC 
COMPLEX OMEGSTR* OM£GCOR. BETASTR* WORWO ♦ SLOPE 
COMMON / NORM / B(3) , TAUO) * NVEC . TAUVEC * PRINT 
COMMON / ROOT / BETA I? . BETASTR ♦ SLOPE 

F (N?*NVEC*TAU?.TAUVEC*SIGN> = N?*N?*(NVEC*NVEC-1.) ♦ TAU?*TA1J?<MTA 
SUVEC*TAUVECM.) ♦ N?*TAU?*NVEC*TAUVEC*2. 

REWIND 5 
READ ( S *501 ) 

REAO<5*500) EPSILON 

B ( 3 ) = B ( 1 > S TAUO) = T AU ( 1 > 

30 OMEGA? = CMPLX ( 0*0 *0.0 ) 

SOLN = 0.0 

BET A I ? = BETASTR 

BETA 12 = BETASTR - OMEGA?*SLOPE 

OMEGAC = OHEGA*EPSILON*EPSILON-»OMEGA? 

BETA IC = BETAl*E°SILON*EPSILON*BETAI? 

ANC = l./GAMMA-BETAIC/AMACH 
CALL NTAU < ANC . OMEGAC ) 

20 CONTINUE 

omeginc = -.10 

OMEGA? = CMPLX ( -10.0 . 0.0 ) S OMEGSTP = OMEGA? 

SOLM = 1.0 

PPINT = 10HNO PRINT 

KOUNT = 0 

IVEC = 1 

21 BETAl? = BETASTR - OMEGA p*SLOPE 
OMEGAC = OMFGA*EPSlLON»FPSILON»OMEGA? 

BETAlC = BETAI+EPSIL0N*EPSIL0N*BETAI2 

ANC = i./gamma-betaic/amach 
CALL NTAU ( ANC ♦ OmfGAC ) 

GO TO ( ?2 . ?3 ) * IVEC 
2? G = F(B-B(2> .NVEC*TAIJ-TAU(2) .TAUVEC. SOLN) 

IVEC = 2 

OMEGA? = 0MEGA2 + CMPLX ( OMEGINC . 0.0 ) 
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GO TO 21 

23 G1 = F (B-B (?) .NVFCt TAU-TAU (?) ♦ T AUVEC • SOLN ) 

IF ( Gl* Gl.LT. l.E-18 .OP. (Gl-G>**2 #LT. 1.0E-10Q > GO TO 24 
6ETAI2 = (OMFGSTR*Gl-OMEGA2*G)/(Gl-G) 

OMEGSTP = 0ME6A2 t OMEGA? - BETA 12 

KOUNT = KOUNT ♦ I 
IF ( KOUNT .GT. GOO > GO TO ?5 
G = Gl 5 GO TO ?1 

24 CONTINUE 

WORW0 = OMEGAC / CMPLX < 1 . 84 1 1 8378. 0 . 0 > 

WRITE (6.606) 

WRITE (6«600> RET4I2, 0ME6A2. EPSILON. OMFGAC. WOPWQ. BET A IC » ANC 
PRINT = 10H % WRITE (6. 60S) B(l) . TAIJ(l) 

niS° = SORT ( (8-B (2) ) * (B-B (2) ) ♦ <TAU-TAU(2> >* (TAU-TAU(?) ) ) 
WRITEC6.608) OISP 
WRITE (6.701) KOUNT 

N? = B - B ( 2) $ TAU2 « T At J - TAU<2> S WRITE (6. 70?> N?. TAU2.G.G1 
IF ( SOLN ) ?B . 20 . 26 
28 EPSILON = EPSILON ♦ 0.10 

IF ( EPSILON .LT. 0.51 ) GO TO 30 
8(1) = 8(3) % TAU(l) = TAU(3) 

RETURN 

25 WRI TE (6.607) EPSILON 
GO TO 24 

26 OMEGA2 = CMPLX ( 10.0 ♦ 0.0 ) 

I VEC = 1 * OMEGINC = -.10 

SOLN - -1.0 * KOUNT = 0 

PRINT = 10HNO PRINT $ GO TO 21 

500 FORMAT (25X.F10.0) 

501 FORMAT!//) 

600 FORMAT (*0 8ETAI2 EQUALS *?G21 . 14. 1 2X* OMEGA? EQUALS *?G21.14// 

S* EPSILON EQUALS *F10.S// 

$* OMEGAC EQUALS *2G21 . 14 ♦ 12X* FREQUENCY RATIO EQUALS *2621.14// 

$* PETAIC EQUALS *2G?1. 14. 12X- N CORRECTED EQUALS *2621.14) 

605 FORMAT (*0 THF INTERACTION INDEX IS *G21.14.?3X* THE SENSITIVE TIME 
S LAG IS *G2 1.14) 

606 FORMAT (*0*135(IHS) ) 

607 FORMAT (*0*SO ( 1H*) * FOR EPSILON EQUAL *F4.2* IT DID NOT CONVERGE*) 

608 FORMAT (*0 NORMAL OlSPLACFMENT = *G21.I4) 

701 FORMAT (*0 KOUNT EQUAL *110) 

702 FORMAT (* N2 = *G21.14* T AU2 = *621.14* G = *621.14* Gl * *621.14) 
END 
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SUBROUTINE WRyTE 

integer tape NO 

COMMON / BLKA / A <3 1 ) 

COMMON / RLKB / B(367) 

COMMON / BL<C / C(^1) 

COMMON / BLOC* / D<IB05> 

DATA TAPE NO / 1 / 

GO TO ( 1 ♦ ? ) , TADENO 

1 T APp NO = 2 
REWIND 1 

WRITE ( 1 ) (A (T) '1=1,31 ) t (8(1) *1=1 *367) ♦ (C C I ) « I =1 ♦ 363) 
$=1.1805) 

REWINO 1 
RETURN 

2 TAPE NO = 1 
REWINO ? 

WRITE (2) (A(I> 'I=l*31> * (8(1) '1=1,367) • (C ( I ) . I =1 . 36 J) 

$= 1 . 1 805 ) 

REWIND 2 
RETURN 
ENTRY REED 

GO TO ( 4 , 3 > , Ta*>ENO 

3 REAp>(2) (A(I) .1 = 1,31) * (8(1) *1 = 1 *367) ' ( C < I > *1 = 1*363) 

$ = 1 * 1 805 ) 

RETURN 

4 READ ( 1 ) (A(I) ,1=1,31) , (B( I) ,1=1,367) * ( C ( I ) » I =1 , 363 ) 

$= 1 » 1 8CS ) 

RETURN 

END 


. (0(1). I 

(0(1) .1 

(DU) .1 
(0(1) .1 
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THE FORM OF THE INPUT DATA IS 


REFERENCE 0.930 0.00 3.700000000.330000001.20 0.9166 

5. 0 0.0 0.0 0.90 0. 010 0.010 1 1 0 

0. 010000000. 000000000.900000000. 00000000 
THE VALUE' OF EPSILQH IS 0.10 
- 20.0 + 0.10 + 10.0 


0.00 
3 3CH0 



THIS IS. ENGINE NUHFtCH REFERENCE 

THE HACH NUMRER 15 . 33000000000000 

THE FREOtlENCT IS 1.7307127532000 0. 

THE PRIMARY M0[>c ASSUMEO IS LHAT a 1 MHAT ■ 1 NHAT a 0 

W/MO a .9*000000000000 0. 

HETAN ■ 2.74999999999995E-02 0. 

PETAC a .16666666466667 0. 

RETArI • ,22841574691221 .12798019395453 

•RETAwlD « >22841574691221 .12748(119395453 

RETAII I) ■ -.12561752604581 5.42142405297 I68E-02 

RETAIC 21 a -B.fifc37l*202?7<»20E-02 .23427026213829 

BETA] I II a -4.740635B9737710E-02 .21136353658143 

8ETAII 41 m -5.61 807 1571 29248F-02 .20474951848918 

RETAU SI » -5.30S123717S814SE-02 .2(508343588394 

RETAU 4) * -S.l 0830488334T57E-62 .21411940711915 

8ETAIC 71 a -5.10049049770763E-02 .21421041174020 

THE INTERACTION INDEX IS .T0T09480293063 

THE INTERACTION INDEX IS <603306438461 13 


THE COEFFICIENT HATMIK IS 3 X 30 
THE RATIO OF SPECIFIC HEATS IS 1.2000000000000 
THE LENGTH OF THE COhHUSTOR IS 2.7000000000000 
THE LINER BEGINS AT 0.0000 ANO ENOS AT .9000 

0 a ,91464664404406 0. 

K a 5.0000000000000 0. 

N a .14116440329633 -.38781876955917 I 

N « .14116440329633 -.38761876955917 I 

N « 1.2139925031691 -.16430072887793 | 

N • 1.1019307334024 -.70990986524754 I 

N ■ .97698896658718 a. 640495559 31767 I 

N ■ 1.0035779264028 -.63560660299751 I 

N a .99409465810853 -.65176798752709 1 

N ■ .48813045101053 -.648B4668623986 I 

N a .98789971205174 -.64912306593949 I 

THE SENSITtVF TIKE LAG IS 2.486635323S4I6 

THE SENSITIVE TIME LAG IS 3.226999074ol39 
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A A A A A A A A A A A A A A AA A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A AAA 
SECOND ORDER. NO LIN^r 

VVVVVVV VVV VV VV W V VVVVVV vvv vvv vvvvvwvvvvvvvvvvvv vvv 

N L EQUAL 1 TO 1 M CQu^t o 

0 5*4765ft6€-0:7 ?i623474E-02-7>176ft64F-92-4.R2?>n5F-03 6 0 9596468-04 4#3]9266F-0 4 

1 2.*0197«E-0? ) .6560?0F-02-2.77O54Sf-O? !.U#l276f-0? 7 0 94264?E-04 ~a , 749603F-05 

2 2.570?9]E-02 1.655774F-03-9.595699F-03 l . 00R094F-02 3 0 I 5972PE-04-? , 56975 ) £ -04 

3 4. *if»42?7E-0?-7«637??4E-fn-3#4007inF-fl3 4,206667F-n7 1 0 1 2797 7E -04-1 .5612R?E-04 

4 -1 .79<.550€-0? 9.224667F-03-] .H7B304F-03 1.64395AF-03 7a ?39?g6F-nft-9 * 77 14 17F-05 

5 -5.1761 | 2E-03 l4)271F^03-9 # SlS7?ir-04 6- A] 20 70F-04 4 a R59P63F-Q5-6 .53RZ49F -oS 

6 -2.4»0®9 1F-03 O.«41346F-<u-s.A59.7|*f- 04 3*)81605F-O4 4 o ?2|376E-nS-4.773??4r-0S 

7 -1 .S«9SOSC-03 6.19S4Q2F-04-7.345497P-04 l*7|404?E-04 3„ ?7 743TC-0S-3 ♦ So 1 692F-05 

fl -1*1365008-03 7.9ft5b?9 F-04-7.3153S2F-04 9.754I64F-Q5 ?, 97177SE-05-2.66 ] 02?F-0S 
9 -5.15S179E-04 7«055441F-Q4-1.47?]37r-04 6.5405208-05 2 . 39736?F-f)5-? . 0 76095F-05 

10 -6.44«9]2E-04 ?.2l4ft5bF.n4-l.l?fl0H4r-o4 4 . 26*1 6?F-05 ? , 1 95692E-05-1 . 6?47 1 6E-0S 

11 -4.9*67228-04 l*BS03O0F-O4-7*741S44r-0S 3.295727F-05 1 . A04652F-05-1 . 3??S9lf-0^ 

12 -4.19770PF-04 1.426167E-04-6.294S27F-05 2.34H24F-05 1 , 6673948-05-] a 06061 l E-OS 

13 -3*7771758-04 1 *?5l 35bF-04-4.697033r-05 1 « 975064C-05 1 . 35A604E-05-5 . 96Q944F-06 

14 -2.967767E-04 1 -»10i497F-O4-3#944O00F-0S U4fll*?7F-05 0 . 29001 2E-0S-7 * 7397A \ £-06 

15 -2o4S7974E-04 9.07O3f,8F*0S-7*0Ob«;oor-o5 1 . 3214768-05 1 . 191*678-05-6. 39675SF-06 

16 -2.21411*8-04 7.446916F*<)5-?*ft716ft9F-05 1.0?ft43ftF-05 1 . i??R44F-05-5 • 33 1 6948-06 

17 -1 .*754768-04 6.496O52F-OS-2.10O5A4F-O5 9.51 34048-06 5 «, 7625S6F-04-4 . 74775SF-06 

IB -1 *7215468-04 5*7&7562F-05-]*922fl4?''-05 7.5971328-06 A*274R438-06-4 - 030S9?f>06 
19 -1.4«04328-04 5.42gbQ9r^ 0 5-W545439F-05 7.20F739F-16 7* 1 634P5F-06-3 . 6flfl | 0 4fc*-06 

?0 “l . 3776R0F-04 4.6054958-05-1.446722^-05 S.566417E-06 6. Rl 4S40E -06- 3 . l 4 75078-04 
?l -1 e 1995458-04 4.39J3048-05-1 .151145F-05 5 • 669 3698-06 5* 95?<*ft38-n6-2 * 92 1342E-06 

22 -1*12*4508-04 7*7662298-05-1 .127574F-0S 4 . 6A 1 ?3 AF-66 5.6999768-06-2.5279648-06 

23 -9.9276618-05 7*6277468-05-9*351904^-06 4.5569978-06 5. 0 1 74A78-16-? . 3774 178-06 

24 -9.4179918-05 1# 1 393A7F-05-9. 0 3750 7F-06 3.R3C95SF-04 4 • A326948-06-? - 06H6678-06 

25 -5. 350G62F-0S 7* 049 1 678-05-7 .57 AS 14F-06 3.794454F-06 4 . 2A254A8-06-1 . 97 1 4298-06 

26 -7.952SM2F-05 2*6543758-05-7.4 1 *0478-06 3 o 1 95 3678-06 4 . 146 1 528-06-1 . 7265908-06 

27 -7.1257978-05 ? *5998458-05-6, 270270«'- O 6 3*1942AAF-06 3. 695469E -06-] .661 61 Af-06 

2fl -6.554? 19E-05 2.250*1*8-05-4. l9ftfW>4F-06 2*713ft56F-06 3, 594 2A3E -16-1 * 46 127 \ E -06 
29 -6* I54fll 68-05 2 .2437408-05-5. 2773Q1F-06 2*730430E-O4 3# 21 9931 E-06-1 #4) 9921 F-06 

N L EQUAL l TO 1 M EQUAL 2 

0 5. 490 7? 7 E -02-5* 5 667 55E- 03 9.745175P-04 5* 565 55 0E* 04- 1 • 75145] E-0 4-1 #?57|)70E“04 

1 4.509200F-02-2.60254RE-02 1 *099791^-03-1 • 052563E-04-2* 1 1957 IE-04 6.509926F-06 

2 1.799Qs7E-0?-l*5h4632F^02 4 , 3359] RP-0^-1* 571 572E-04-5 . 524 ft5«F-n5 6.7?Hfl5?F-05 

3 6. 11 513 1F-91-3* 1575335-03 1 • S39769P-05-2# 1394 33F-04-3 . 1 5 71 17F-0S 4.406 729E-05 

4 1.575443F-03 5*8975?9F-04 9.fl94?16r-o5-l , 3) 9623F-04-2 . 0 707S7E-05 2 * 9757H5E-05 

5 -1 • 367031 E-05 9.4A2H3F-04 6 . 6 t P«b9r-o5-R • 659 75 7F*05- U 3991 94E-05 2.166552F-05 

6 -4.SI4522E-04 4.412530F-04 5. 691546F-05-6. I P4957E-05- 1 . 2666Q9E- 05 l#7?S14flF-0S 

7 -5-I003J5E-04 4.6359?0F^ft4 4.374] 77F-05-4 *5] 0752F-05-1 . 0?l470K-95 U39JH03F-05 

6 -4.9<*3sp2£-04 -U 1607HF-04 3 . 926 1 62f-rS-3 • 379 1 35F- 05-9 . 933950F-06 l*16?5?9E-05 

9 -4,?6??2?£-04 ?*5776l7E-04 3. I359**OF-05-2.6U 54 1F-05-5.44S9 llE-06 9.712204F-06 

10 -3.«*2999?F-04 1.9652^0F-04 ?. 54^551 F-05-7. 023179F-0S-5. 291 ) 73F-06 5.237n?F-06 

11 -3.233293E-04 1.609349F-04 2. 3?n9rt3r-o5-l *639605F-0S-7 , 1 6931 3F-04 7.029550E-06 

12 -2.5946?2F-0<* J.250326E-04 ? • 1 24095^-05- 1 • 30626 \ F-05-6 . 997217E-06 5.952539E-06 

1 3 -2.463993E-04 1.09681HF-04 1 . 764474F-05-1 . 1 0 ] 7Q4f>05-6.0779?0E-f\6 5.155145F-06 

14 -2.?319k 5F-04 5. 9H 350 7F-05 l .671 031 F-oS-5.9651 50F-06-5 .91 0 1 1 0F-06 4. 435717^*06 

15 -1 . 91 5979E-04 7.95746 1F-05 1 . 37577 1 f- 05-7. 53] 544F-06-5. 1 54M56F-n6 3#910R73f-n6 

16 -1.755655E-04 6*656349F*rtS 1 . ?P4?5QF-0S-6« 5ft 7994F-Q6-4 . 9Q97 0QE-06 3i3660l6F-n6 

)7 -1.525916E-04 6*«429l4F-05 1 * ft97570F-05-5. 6246S5E-06-4 . 3794t9£-06 3.012496E-06 

16 -l*416f)70E-04 5 » I 75376F-05 1 . 033672F-05-4.913579E-06-4.2<»4477e-04 2.6094R0E-06 

19 -l.?4l1??F-04 4. 749ft 5 1C— 05 ft . 93431 6F-06-4. 49235RF.-06-3 . 73794 1 E-06 2.37Q645F-06 

20 -1.I42IA2E-04 4.066240F-05 5 . 45 756] F-Ob-3 . 835M RE-06-3. 625442E-06 2.066169F-06 
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21 “1 .429195E-04 3.H34949E-05 T.400S40f-06-3.56762TE*06-3.?UV434t-06 I • vo^vjlt-04 

22 -9.4969T9E-0S 3.331521E-05 7. 0766STE-06-3.0T60T3E-06-3. 1 2169AE- Oh 1.646162E-04 

23 -6.450O2SE-0S 3.163234E-05 6.227?30r-06-2.90l69?E”06-2.774449E-06 1.5S5401E-06 

24 -0.2046R9F.-OS 2.7699O7E-05 5.9077Hir-06-2.S220l6E-06-2.705?60E-06 1 .370469E-06 

25 -7.366036E-05 2.4S51S5F-05 5.30l903 r -06-2.40703?E-06-?,415310E-06 1.292O72E-06 

26 -T.91170Se-05 2.3405J.3E-05 S. 1 2««»47F-06*2. 10*«9OOE-0R-2. 341 34SE-0R 1.147943E-06 

27 -4.34T760F.-05 2.261294E-05 4 .S60542E -04-2 . 029491 F-06-? . 11 72C6E-06 1.Q686UE-06 

28 -6.4R9029E-05 2.0046C3P-05 4.441S7TF-06-1 .TR55TSE-06-2.075419E-06 9.679765E-07 

29 -S.52373IE-05 ).96961lf>05 3.975921 E-06-1 .734942E-06-1 *848 166E-06 9.291971E-0T 


THE SLOPE 1$ 26334236*21269 1.3660352949050 

THE tNTFROCTION INDEX IS .6419961108*035 i THE SENSITIVE TINE LAG IS 3.2566 2 2051 3033 

imtmimmfiTiiimmmfitmmtumiimmfmmmmiiMiittmmmmiiimttftmttmimmmmitmimmmmi 

PET A 12 EQUALS .55576107101676 -.74430425472799 0HE0A2 EQUALS .492029147-9944 0. 

EPSILON EQUALS .10000 

I 

ONEGAC EQUALS I . 7356 130444750 0. FREQUENCY RATIO EQUALS .9424723521727* 0. 

6CTAIC EQUALS .23397335742239 .12031713140725 N CORRECTED EQUALS .12432315972006 -.36459736790074 

THE INTERACTION INDEX IS .59676136399513 7HE SENSITIVE TIME LAG IS 3.2414331043256 

NORMAL DISPLACEMENT « 1 .564036S4279TT6F-02 

MOUNT EQUAL 22 

N2 ■ -6 ♦ 525 0 74465994 15E- 03 TAIJ2 * 1 .443402971 1 7651E-02 E. « -1.S2631505220372E-09 Gl • -5.B37653n922fl429E-10 

RETA12 EQUALS .55596459735059 -.76694796094676 OHEGA2 EQUALS .49249621669194 .0. 

EPSILON EQUALS .10000 

QMEC.AC EQUALS 1.7356377353069 0. FREQUENCY RATIO EQUALS .94247489963477 0. 

PETAIC EQUALS .23397459206572 .12031071434504 N CORRECTED EQUALS .12411941549763 -.36457792225769 

7HE INTERACTION INDEX IS .59673656200702 THE SENSITIVE TIME LAG IS 3.2414161702469 

NORMAL DISPLACEMENT * 1 .S6426S09776367F-02 

•COUNT EQUAL 26 

N? * -6.5670764541 o?6 2E-0 3 7 AU2 * 1 .44170934350609E-02 G « -1 . 2*575602666350E-09 Gl » -4.75636561S59390E-10 

AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

SECOND ORDER* R/ LINFR 

v V V V vv V V V V V WV V V V V V V vv V V V V vvv V v v vv V V VV V V V WV VV V VV V 

N L POUAl 1 TO 3 H FOUAL 0 

0 5*IO<»S4?E-01 ?• 7<**Z>76E-fl2-W41 5* 745*<n£-04 1 . 095A31 E-03 

1 2.7A<M?1E>0? ?• lA2f>A3F-f)?-? # <M?37)>r-02 Ai445*«*H-03 f>*?74me-0<i W07470AC-03 

2 2*?1?<»0SF-0? S*4A72fWF-<n-1.0SQ7A2F-O2 B,447747f>03 3« \ \ 23E-04 4 • 371469004 

3 7.bSl?23«t-<)l-4. A9440?^^03 3# 5*A3R0F-A3 ?• A04ft*#»E-0<* 1 . O30347F -04 

4 - \ .4linS0£-0?-S.03;»?Q o r-r> w # *4?!SVF-n3 1 .4RA4*4F-0;* 3*0?4fll7F-04-7#l 1 7A1SF-05 

5 -3.1?7?07£-0 7-?. !494^7F-rt3-i # pQiS2Mr-0 3 **.A<SVi1«iF-04 ? .09561 «F-ft4-A*H9 111 ^-05 


C-61 



e ~? B n39oiAe°04“4°44969?F-o4-i a?i5?j4F~04-6ortP30fl4E°06°9o9;i«>im-04 4jo<m^<H.°os 
_S -4*.4A424lE«*04-4*9?7i»71F»O4-l ^SftSAF-QiW^OAayfclE-Afc ? * 76009f,£ <*AS ] 0 954K?oF-05 

10 -3. 483?67£-0*-4.057S07F-0*- 1.421720^-04-8. 2S7O9SF-06 5*491 946F-05 5, 67 3342F-06 

1 1 -2.*1A?59E-O4-3.2?O?71F-04-l.On?S63F-O4-6.79ft5ft5F-O4 4.431451E-05 4 . ?94080F-n4 

12 -2.3O46O2e-O4-?.769272F-O4-A.n?3ft9ftF-O5-l,O5O04?F,-GS 1 ** 7 9,?B?497F-06 

13 -1 .8*71 32E-04-?o273979E-04-2.6Q??09r-04-l .OlASS I F.-05-4 . 1 5 1 f 0 l E-06 1 .22H 12QE-0S 

14 -1.465?l?e-O4-?*d05ll ]F-04-2*277739F-05-1 . 1 3747 OF-O 5-3* 321151F-06 1 .163934E-05 
J5 -W783664E-04-1 .6BA454F-04-?.5943A]F-O5-9* 191915F-06 R % ?4 1 8A9E-06 *,44\l77F-06 

16 -i.)55054F-04-| .S2O646F-04-1.1594O1F-O5-8.7676ROF-O4 ) .7S7747F-0S 6,11P100F-06 

17 -9.193182E-0S-I * 10491 3F-04-2.50167OF-05-6. 92*0*0^-06 1 .496740F-05 5* 0<>07H?E-06 
IB -9.16A937E-0S-1 . 197SS4F-fl4-1.654994r-05-6.973356F-06 6. ?9?5 1 8 F>o* 5*451 IH3E-06 

19 -7* A| 7907E-05-J .O41037F-O4-7. 7tn395F-06-S.4Rf)734F-06-l *042?A4E-06 5,544 169 E -06 

20 -7.493954E-OS-4*656RlBF-OS-7.549912F-06-6.16031^F-06-6.A30S74f:-n7 5. 34S>?5SE-06 

21 -6.144959E-05-A.497411F-05-9.37?S49P-06-5*236974E-06 3.79768SF-06 4.486923F-06 

22 -S.A3?579E-0S-7.9S07S3 r -OS-l*?32062F-05-5.149A0IF-d6 7.79I919F-06 3,«7H?39F-06 

?3 -4*94?4S2F-05-7* ^7234 ir-<)5- 1.024 U A^-05-4*2H06 7?F-04 6.A39449E-06 1. 1680R?F-06 

24 -4.9605?OE-0S-4*475075F-0S-7.?44] 2SC-06-4. 247462F-06 1*?11449f-04 3.329133F-G6 

25 -4.339] 77E-05-S.980756E-05-3. 64? 79?r-06-3.66?267E-06-l .177443F-07 3.17492M-06 

26 -4.?Hfl4lS£-05-5*6767l £, F-05-l* 7BOOBAF-06-3.744 79AF-06 5. 99*4 l frE-nft 3.09*9}4F-04 

27 -3.42063SE-05-5. 1 i486 IF- 05-4. 697 67 2F- 06- 3 .267800 E-06 ?. ] 78681 E-06 2.751 49?F -06 

28 -3.534?69E-OS-4.H84776F-05-6.3?49?2F.06-3.2752Q*F-06 4*2173?ft£-06 ? .542447F-06- 

29 -3. 06] 065E-05-4 . *302635-05-5. 3RG067F-06-2.803719E-06 3.T464A1E-06 2.264446E-06 

N 

a 
1 
2 

3 

4 

5 

6 

7 

8 
9 

io 
u 
12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 
73 

24 

25 

26 
27 
?fl 
29 

NVEC * -.49847661286642 TAUVEC » 5 . 60740467604330E-02 


THE SLOPE TS -.50*S0l4258852a 1.4*20197171086 

THE INTERACTION INDEX FS *70583873874468 THE SENSITIVE TIME LAG IS 2.4844325867333 


L EOUAL 1 TO 3 m EQUAL ? 

5.27251 SE-02-5. 95654 1F-03 1 . 14447 OF-03- U298522E-04-2 . 863508F- -)4 2 o273800F -04 
4.)93529E-02-2*Z06012F-0? 1 .3S723OF-03-1 .31 651 *£-0 3-3.8071 07E-04 4 * 1 2?5R7F-(I4 
1.517787E-02-1 .58545*5-62 4.4771 8Sr-64-l , 0*44 7 Af -*3-] .307637F-04 4o29?6S|F-04 
4. 79235 3E-0 3-4*53888 1 F-0 3 1 . 4?62 iSF-QS-3 .937355F-04 4.724?Q5F-n5 lo576fll?£-04 
1.9354*2E-03-2.8fl2926F-n3-l.2373K2F-04-4.925419F-05 1 . 06758 1 f-f)4 6 o 370954E-07 
6.449676E-04-|.548242F-03-8.426Apsr-D5 2.1 75480F-O5 7, 7 1 6002E-05-? . RS3547E-0S 
-l.n*,374ie-0*- t >*43l4*,9£-04 3 . 2 7 720 7 F - 05“ 1 • 86 0 I 8 IF - 05 4,29?H70C-06 6 o 683?0SE-06 
-3.179S01E-04-6*4fi0821F«04 9. 939204F-05-4 .550 79 2F-05-4 . 2S45S3E-0S 3.531 774E-05 
-3.224479E-04-9.9H6870F-04 8 .60B664F-05-3. 279957E-05-3, 967 1 02E-05 3*11 74UF-05 
-1.914437E-04-4.004164F-04 ?*S5342?r-05-S . 90BO8 7E-06-4 ♦ 79] 433E-06 R.607683E-06 
-1 . 2*£i5nOE-04-3 . 186 O 6 OF- 94 - l*?l9840 r -05 1 • 119999F-Q5 ? . 0 1 88O7E-05- 7 . 9] 2454E-06 
-’l.^69778F-04-?*B08008F-04-1 . 1537^-05* U?74|93f-0S 1 .76] ) \ 5E -05-8* 99541 5£-06 
-l • 36 3841 E-04-2*40S09 I E-04 1 .42575RF-0S 7. 81 9O40f>O6-3*94OO7SF.«O6-| c512B99F-06 
-1 .14991 9E-04-?.01?4R4F-04 3 . 0 0 0 0 S5f- 0 5 *3 1 4*9* 1 8 3F-06- | ,93?4B6f-05 4. 1R3278F-06 
-I.268804E-04-1 .T63056F-04 ?.665R4?r-05 3 .5S2444F -06-1 . 744943F -05 3 .6] 9350^-06 
-9.B40700E-OS-1 .513143F-0* 4. 63^57^66 5l245'ri5F-06-3 *8389 77E-06-3. 90941 4FT-07 
-7.3?9397E-05-1.36<.198F^04-1.765013F-06 6.937439F-0* 6 . 35070HE -06-3. S50601F-06 
-6. I45008E-05-1.19A244F-04-2.045170F-06 6. 643*>V3£-04 5. A63993F>06-3 .631 57 3E - 06 
-6.87534RE-05-1.089392E-04 6.376049F-O6 5. 738332F-06-2. 6 1 04B8F-06-1 .484? I 7F-06 
-6.41 3983E-05-Q. 5944| 3F-05 1 . 194867F^3 4 . 499974E-06-8. A83871F-06-S* l 662S4F-07 
-8.338706F-05-B.81AOH5F-05 1*1081 I3 r -0S 4. \9R249£-06-A.24?394E-ft 4-4. 97239] Pf>07 
-4.9aS569F-05-7#8.1942lE-05 *.621 986^-0# *4 . 069 1 74E-06-? . 1 ft 7ft AQ£-n6-l .279467F-06 
-4* 29Q264E-05-7* 3033 5ft F-05 9.3?Afi5?r-0ft 4.285704F-O4 2 *4?53? 1 £-06-? • 0 38442F-06 
-3.691 942E-0S-4.571644F-05- 1*87474 IF-bT^ 4. 00O?5*F-ft6 ? . 3359S6E-06-? , 0005A0E-06 
-4. 02721 Of -05-4. 1 73597F-0S 3.50788 7F-06 3 . 748 0 1 7E-04-1 . 61 4H64E-06-1 .544S92E-06 
-3.831 508E-05-5 *5 7882BE-05 6.070357F-04 3 .Z3S554f-06-4 , 63 1 7 1 1 E-04-] o 034974C-06 
-3.763743E-05-5.259042F-05 5.7*«l«?F-06 3.06424 Af-06-4 . 381 4) 2F-06-9, 777572F-07 
-3.«34194E-05-4.773l«6E-05 7.6755S3F-06 2.8 3323»E-04-l • 333939E-06-1 * 1 23779E-06 
-2.772258E-05-4.534255F-05 4.8H302F-07 2. 8462 1 7F-04 ] . 045595E-06- 1 * 34546 AE-06 
-2.428152E-05-4.148637E-05 ?.B]3373F-07 2.64B50BE-06 1 • 04B75 IE-06-1 • 29S646E-06 
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\\*\ ** r> ** fr*% Vt*> *'<>*> l » 


QETAI2 EQUALS -5.693979277721A8E-02-6.7I601683625633E-03 
CPS t LON EQUALS .10000 

OMEGAC EQUALS I .7302006921376 0. 

METAIC EQUALS -5. 157630290*Rfc«*E-02 .21*1 *3*515918* 

THE INTERACTION INDEX IS ,707560176731*9 
NORMAL DISPLACEMENT ■ 4.73*1 I246O3065SE-0* 

MOUNT EQUAL 20 

N2 » 4.7337380086IT6SE-04 T*U2 - -5.95*?659R69?939E-06 


OMEGA 2 EQUALS -5. 12061062* 3R897C-02 n. 

ERFOUENCy RATIO EQUALS .939721 flft*8721 0 0. 

N CORRECTED EOIIALS .98*62516031772 -.6**9) 995027*29 

THE SENSITIVE TIME LAG IS 2. *666293693257 

» “*.2* 322386 1 906* IF” 1 0 


0 « -1.11089916A20696E-09 61 




METAI2 EQUALS -5.6**102*2*T3706E-02-0.181 19563079991E-03 

EPSILON EQUALS .10000 

0ME6AC FOUALS I .7302105785027 0. 

RETA1C FOUALS -5. lS713lS2l95*9*E-02 .21*12879900380 

THE INTERACTION INDEX IS .70753*75591*33 
NORMAL DISPLACEMENT ■ *.*2600 1670 l*3SHF-0* 

MOUNT EQUAL 21 

N? » 4..399529R3700920E-0* TAU2 ■ -*.»335085651?7*3E-0S O ■ 


OMFOA? EQUALS -5. 02 1 7*6977?08T0F -02 0. 

FREQUENCY RATIO EQUALS .93972725***207 0. 

N CORRECTED EOUALS .9096100*61196* - .6*8075 15097060 

THE SENSITIVE TIME LAO IS 2.4B65B6988S060 

-1 .*56779*1 3420B2E-09 pi ■ -5.S64*S2*5766*8SE-10 





